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Abstract 

-^^efect  properties  of  copper  are  calculated  using  molecular 
statics  with  an  interatomic  potential  recently  derived  from  first 
principles.  Tri-  and  tetravacancies  are  found  to  be  very  mobile 
with  migration  energies  of  0.56  and  0.39  eV,  respectively, 
compared  to  previously  calculated  single  and  divacancy  migration 
energies  of  0.82  and  0.55  eV,  respectively.  Using  the  binding 
and  migration  energies  calculated  with  the  interatomic  potential, 
annealing  kinetics  in  copper  and  modeled  using  rate  equations. 

The  effective  activation  energy  of  annealing  in  the  model  is 
within  0.02  eV  of  single  vacancy  migration  energy  over  a  wide 
range  of  sink  and  initial  single  vacancy  concentrations,  which 
conforms  to  experimental  results.  In  two  cases,  however,  the 
larger  clusters  affect  the  activation  energy  and  no  definitive 
conclusions  about  whether  or  not  the  calculated  cluster  migration 
energies  are  correct  for  copper  can  be  made. 

The  stability  and  structure  of  larger  vacancy  clusters  with 
ten  to  forty  vacancies  were  also  investigated  using  the  first 
principles  copper  potential.  The  stacking  fault  energy  was  first 
calculated  and,  for  the  pot^tial  cutoff  radius  used  in  the 
defect  calculations,  yield^  of  a  value  of  65  mJ/m*'  compared  to 
the  experimental  value  of  ^wO  mJ/m* .  To  investigate  the  large 
clusters,  vacancy  platel^s  of  various  sizes  were  created  in  a 
close-packed  plane  and^he  system  was  relaxed  to  the  minimum 
energy  configuration.  QBmall  vacancy  platelets  with  as  few  as  ten 
vacancies  collapsed  into  stacking  fault  tetrahedra  and  faulted 
loops,  depending  on  the  shape  of  the  platelet.  Stacking  fault 
tetrahedra  are  found  to  be  the  roost  stable  large  clusters. 
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Abstract 


Defect  properties  of  copper  are  calculated  using 
molecular  statics  with  an  interatomic  potential  recently 
derived  from  first  principles.  Tri-  and  tetravacancies 
are  found  to  be  very  mobile  with  migration  energies  of 

0.56  and  0.39  ev,  respectively,  compared  to  previously 

calculated  single  and  divacancy  migration  energies  of 
-  0.82  and  0.55  ev,  respectively.  Using  the  binding  and 
migration  energies  calculated  with  the  interatomic 

potential,  annealing  kinetics  in  copper  are  modeled  using 

rate  equations.  The  effective  activation  energy  of 
annealing  in  the  model  is  within  0.02  eV  of  single 
vacancy  migration  energy  over  a  wide  range  of  sink  and 
initial  single  vacancy  concentrations,  which  conforms  to 
experimental  results.  In  two  cases,  however,  the  larger 
clusters  affect  the  activation  energy  and  no  definitive 
conclusions  about  whether  or  not  the  calculated  cluster 
migration  energies  are  correct  for  copper  can  be  made. 

The  stability  and  structure  of  larger  vacancy  clusters 
with  ten  to  forty  vacancies  were  also  investigated  using 
the  first  principles  copper  potential.  The  stacking 
fault  energy  was  first  calculated  and,  for  the  potential 
cutoff  radius  used  in  the  defect  calculations,  yielded  of 
a  value  of  65  mj/m*  compared  to  the  experimental  value  of 
^70  mJ/m*.  To  investigate  the  large  clusters,  vacancy 
platelets  of  various  sizes  were  created  in  a  close-packed 
plane  and  the  system  was  relaxed  to  the  minimum  energy 
configuration.  Small  vacancy  platelets  with  as  few  as 
ten  vacancies  collapsed  into  stacking  fault  tetrahedra 


and  faulted  loops,  depending  on  the  shape  of  the 
platelet.  stacking  fault  tetrahedra  are  found  to  be  the 
most  stable  large  clusters. 

In  order  to  investigate  the  range  applicability  of  the 
potential,  molecular  dynamics  calculations  of 
thermodynamics  and  vibrational  properties  were  performed. 
Parallel  calculations  were  also  done  using  an  empirical 
potential.  It  was  found  that  the  agreement  between 
calculated  and  experimental  properties  was  much  better  at 
low  temperatures  than  at  temperatures  near  melting  for 
both  potentials.  Neither  potential  was  found  to  be 
significantly  better  than  the  other  for  thermodynamic 
calculations. 

A  first  principles  potential  for  silver  was  also  used 
to  calculate  defect  migration  and  thermodynamic 
properties.  The  divacancy  is  found  to  be  the  most  mobile 
vacancy  cluster,  followed  by  tetra-,  tri-,  and  single 
vacancies.  Combined  with  the  copper  results,  it  is 
concluded  that  di-,  tri-,  and  tetravacancies  are  all 
highly  mobile  compared  to  single  vacancy,  and  the 
relative  order  among  the  clusters  depends  on  the 
potential.  As  with  the  copper  potentials,  the 
thermodynamic  properties  calculated  using  the  silver 
potential  were  in  good  agreement  with  experimental 
results  at  low  temperatures,  but  discrepancies  were  found 
at  high  temperatures. 
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Introduction 


The  technological  importance  of  vacancies  and  small 
vacancy  clusters  in  materials  has  led  to  a  considerable 
research  effort  toward  studying  their  effects.  For  many 
reasons  this  has  not  been  an  easy  task.  The  primary 
difficulty  is  the  small  size  of  these  defects  which  is  on 
the  order  of  a  few  lattice  constants.  This  makes  it 
difficult  to  identify  and  distinguish  from  among  vacancies 
and  small  vacancy  clusters  and  has  often  forced 
experimentalists  to  use  macroscopic  properties  to  quantify 
defect  concentrations.  For  example,  single  vacancy 
migration  energies  have  been  estimated  for  most  metals 
using  annealing  experiments  in  which  the  electrical 

resistivity  is  used  as  a  measure  of  vacancy  concentration. 

4 

As  Balluffi  summarizes,  this  method  has  many  drawbacks, 
and  the  single  vacancy  migration  energies  of  some  metals 
have  only  recently  been  reliably  determined  after  twenty 
years  of  experiments. 

In  the  case  of  copper,  there  is  very  little 
experimental  information  about  small  vacancy  clusters. 
Although  the  experimental  values  of  single  vacancy 
formation  and  migration  energy  are  fairly  well 


established  ,  the  binding  and  migration  energies  of  larger 

clusters  are  unknown.  Annealing  studies  done  by  Wienhold 
5 

et  £l  have  suggested  that  at  least  one  vacancy  cluster  is 

as  mobile  as  the  single  vacancy  in  copper.  This  is 

confirmed  by  perturbed  angular  correlation  (PAC)  studies 

6 

of  annealing  in  copper  done  by  Wichert  et  al  and  Deicher 
7 

et  al  in  which  there  appears  to  be  a  second  mobile 

vacancy-type  defect  in  addition  to  the  single  vacancy. 

Neither  of  these  methods  has  been  able  to  specifically 

identify  the  size  of  the  cluster,  although  it  is  assumed 

to  be  a  divacancy  since  theoretically  divacancies  are 

4 

believed  to  be  more  mobile  than  single  vacancies  . 

Defects  involving  larger  vacancy  clusters, 

dislocation  loops  and  stacking  fault  tetrahedra  (SFT),  can 

8 

be  identified  using  electron  microscopy  and  the  size 

distributions  of  loops  can  be  determined  using  diffuse 

9  9 

X-ray  scattering  .  Recently  Larson  and  Young  analyzed 

the  size  distribution  of  dislocation  loops  in  copper  using 

diffuse  X-ray  scattering  and  found  loops  which  were 

apparently  as  small  as  about  10-15A  in  radius.  This 

conflicted  with  earlier  electron  microscopy  studies  which 

10 

found  loops  as  small  as  about  25A  .  Larson  pointed  out 

that  diffuse  X-ray  scattering  should  be  better  for 
analyzing  small  loops,  although  it  does  not  allow  one  to 
image  defects  but  only  produces  a  size  distribution.  The 


diffuse  X-ray  scattering  data  in  copper  is  interpreted 

9,  11 

using  anisotropic  elasticity  theory  .  Although  these 

results  are  no  doubt  more  reliable  than  the  previous 

studies  using  isotropic  theory,  it  is  well-known  that  the 

assiimptions  behind  elasticity  theory  break  down  near 

12 

dislocation  cores 

With  these  uncertainties  in  basic  defect  properties 

of  copper,  it  is  evident  that  computational  and 

theoretical  methods  could  contribute  important  information 

to  the  current  understanding.  The  purpose  of  this  work  is 

to  resolve  some  of  these  questions  about  the  properties  of 

small  vacancy  defects  in  copper  using  the  computational 

methods  of  molecular  statics  (MS)  and  molecular  dynamics 

(MO) .  In  these  methods  the  material  is  represented  by  a 

collection  of  atoms  which  interact  with  one  another 

through  a  given  interatomic  potential.  This  work  is  a 

particularly  timely  application  of  MS  and  MD,  since  an 

interatomic  potential  recently  derived  by  Dagens  using  the 

13,  14,  15 

pseudopotential  method  is  employed 

16 

The  pseudopotential  method  is  used  to  derive 
potentials  from  first  principles  and  with  a  minimum  amount 
of  experimental  information.  Even  with  this 

justification,  however,  the  potential  must  be  able  to 
duplicate  existing  experimental  results  to  be  considered 
reliable.  The  quantities  so  far  calculated  using  the 
potential  derived  by  Dagens,  which  include  the  most  stable 


lattice  type  ,  the  elastic  constants  ,  phonon 
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frequencies  ,  and  the  properties  of  single  vacancies  , 

all  agree  veil  with  experimental  values  for  copper.  In 

this  work  the  verification  of  the  potential  is  further 

extended  by  the  calculation  of  basic  thermodynamic 

properties  and  vibrational  properties.  The  bulk  modulus 

and  phonon  frequencies  are  also  calculated  to  insure  that 

the  computer  implementation  is  correct  since  these  can  be 

compared  with  previous  calculations.  It  should  be  pointed 

out  that  these  results  are  for  nonzero  temperatures  using 

molecular  dynamics  while  the  previous  calculations  where 

done  using  techniques  at  zero  or  low  temperatures. 

The  defect  calculations  in  this  work  for  the  most 

part  follow  standard  methods  in  this  area.  The  first 

attempt  to  calculate  the  displacement  field  around  a 

single  vacancy  was  made  by  Girifalco  and  Streetman  in 
18 

1958  .  This  work  applies  to  BCC  metals;  in  1960 

Girifalco  and  Weizer  made  a  more  complete  study  of  vacancy 

19 

relaxation  in  cubic  metals  using  a  Morse  potential 

Tewordt  used  two  Born-Mayer  potentials  to  calculate 

20 

properties  of  various  point  defects  in  copper  in  1958 
This  work  is  significant  because  it  was  the  first  case 
where  an  atomistic  or  discrete  method  was  combined  with 
elastic  theory.  In  1965  Johnson  reported  on  the  first 
calculations  of  both  the  single  and  divacancy  migration 


21 

energies  using  a  copper  potential  .  Using  a  potential 

for  nickel  (the  Johnson  I  potential),  Johnson  later 

reported  calculations  on  the  formation  and  migration 

energies  of  larger  clusters  and  discovered  a  simple 

22 

relationship  describing  these  values  .  In  1966  Cotterill 
12 

and  Doyama  published  a  work  on  the  investigation  of  the 

edge  dislocations  in  copper  using  a  Morse  potential;  they 

also  reported  calculations  for  the  stacking  fault  energy 

using  Born-Mayer  and  Morse  potentials.  In  1975  Bennett 

discussed  the  method  used  to  calculate  the  migration 

23 

energy  in  static  systems  in  this  work 

The  first  application  of  molecular  dynamics  to 

thermodynamic  property  calculations  was  done  by  Alder  in 
24 

1959  .  Rahman  was  the  first  to  use  continuous  potentials 

25 

for  these  calculations  in  1964  .  The  original 

calculations  were  for  liquids;  Dickey  and  Paskins 

26 

calculated  the  properties  of  solids  in  1969  .  These 

studies  were  performed  on  microcanonical,  i.e.,  constant 

volume,  systems.  Andersen  suggested  a  way  of  doing 

constant  pressure  simulation  in  1980  in  which  the  system 

27 

volume  was  allowed  to  expand  and  contract  .  Parrinello 

and  Rahman  subsequently  did  calculations  using  a  technique 

in  which  the  boundaries  change  shape  as  well  as  volume 

28 

which  permitted  the  lattice  structure  to  change 

From  this  review  it  is  apparent  that  MS  and  MD,  as 


applied  to  defect  and  thermodynamic  property  calculations 
in  this  work,  are  well-established  methods.  Except  for 
the  flexible  boundary  method,  the  general  techniques  in 
these  areas  have  not  changed  significantly.  The  molecular 
statics  method  has  become  computationally  faster  and  new 
methods  have  been  presented  for  calculating  properties 
such  as  elastic  constants.  It  is  interesting  to  note, 
however,  that  there  has  not  been  much  cross-fertilization 
between  workers  in  defect  and  thermodynamic  property 
calculations.  Apparently  the  only  potential  which  has 
been  extensively  investigated  by  workers  in  both  areas  is 
the  Lennard-Jones  6-12  potential  for  argon  and  other  noble 
gases.  This  potential  is  not  as  applicable  to  metals  as 
other  potentials,  however,  and  the  defect  properties  of 
metals  are  of  much  greater  technological  interest  than 
those  of  argon.  In  this  work  a  single  potential  is  used 
to  represent  copper.  The  hypothesis  that  a  single 
interatomic  potential  can  be  used  to  successfully 
calculate  many  different  properties  of  a  material  is  being 
tested. 

The  outline  of  this  work  is  as  follows.  In  Chapter  2 
the  potential  used  in  this  work  is  described  in  detail  and 
compared  to  some  previous  potentials  used  for  copper.  The 
difference  between  empirical  and  first  principles 
potentials  are  emphasized,  and  some  of  the  major 
disadvantages  and  necessary  assumptions  related  to  the 
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potential  used  in  this  work  are  pointed  out.  The 

computational  methods  used  in  this  work,  molecular 

dynamics  and  molecular  statics,  are  discussed  in  Chapters 

3  and  4,  respectively.  The  molecular  dynamics  method  as 

used  in  this  work  is  covered,  including  the  flexible 

28 

boundary  technique  of  Parrinello  and  Rahman  .  A  new 
molecular  statics  method  which  is  an  order  of  magnitude  or 
more  faster  than  previous  methods  is  presented  in  Chapter 

4. 

The  method  and  results  of  thermodynamics  calculations 

15,  1 

are  presented  in  Chapter  5.  Dagens'  potential  , 

derived  from  pseudopotential  theory,  and  the  empirical 
29 

Morse  potential  ,  both  reproduce  copper  properties  at  low 
temperatures  better  than  at  high  temperatures.  The  reason 
for  the  good  performance  at  low  temperatures  is  that  both 
potentials  (especially  Dagens'  potential)  are  particularly 
good  at  reproducing  the  phonon  spectra. 

Chapter  6  is  devoted  to  the  evaluation  of  migration 
energies  of  tri-  and  tetravacancies.  Molecular  statics 
results  show  that  in  the  order  of  increasing  mobility, 
small  vacancy  defects  in  copper  are  single  vacancies, 
trivacanciec,  divacancies,  and  tetravacancies.  To 

investigate  the  experimental  consequences  of  the  migration 
energies,  rate  equations  describing  vacancy  defects  in 
copper  are  used  to  model  annealing  kinetics.  In  most 
cases,  it  is  found  that  the  migration  energies  of  the 


divacancies  and  larger  clusters  do  not  affect  the 

effective  activation  energy,  which  agrees  with 

4,  5 

experimental  observations  .  It  is  also  concluded  that 

the  migration  energies  of  defects  are  sensitive  to  the 

potential  and  are  not  readily  predictable  based  on  the 

results  for  other  potentials. 

Stacking  fault  and  large  vacancy  cluster  calculations 

are  presented  in  Chapter  7.  The  stacking  fault  energy  is 

found  to  have  the  same  order  of  magnitude  as  experimental 

predictions  if  the  potential  is  not  extended  beyond  twelve 

nearest  neighbors.  It  appears  that  the  analytic  form  of 

the  interatomic  potential  which  is  fitted  to 

1 

pseudopotential  results  is  inaccurate  at  large  radii. 

Vacancy  platelets  with  as  few  as  ten  vacancies  are  found 

to  collapse  to  form  stable  and  metastable  clusters.  The 

most  stable  vacancy  clusters  are  those  which  arise  from 

vacancy  platelets  in  the  shape  of  triangles,  rhomboids, 

and  hexagons,  in  order  of  decreasing  stability.  The 

triangular  platelets  invariably  form  stacking  fault 

tetrahedra  (SFT),  the  hexagons  form  Frank  loops,  and  the 

rhomboids  form  a  complex  intermediate  configuration.  The 

displacement  field  immediately  above  the  hexagonal  loops 

is  found  to  differ  greatly  from  anisotropic  elasticity 
30 

theory  ,  which  has  been  used  to  interpret  diffuse  X-ray 
9,  11 

scattering  data 

The  general  pseudopotential  theory  developed  by 


Interatomic  Potentials  for  Copper 

2.1  Introduction 

Molecular  dynamics  (MD)  and  molecular  statics  (MS) 
can  be  used  to  simulate  a  vide  range  of  material 
properties  and  behavior.  Both  micro-  and  macroscopic 
properties  can  be  calculated  such  as  the  vacancy  formation 
energy  and  the  elastic  constants.  Regardless  of  the 
property  being  simulated,  however,  MD  and  MS  both  require 
that  the  interatomic  potential  be  given.  The  power  of 
these  methods  is  that  if  one  knows  how  atoms  interact  on 
an  atomistic  level,  many  other  properties  can  be 
calculated  without  further  assumptions.  In  practice, 
limitations  on  computational  resources  make  it  necessary 
to  restrict  interactions  between  atoms  to  relatively  short 
range,  and  the  the  system  size  a  few  hundred  or  thousand 
atoms.  It  is  well  recognized  that  one's  ability  to 
predict  the  behavior  of  a  real  material  depends  on  the 
interatomic  potential  which  is  used,  and  that  the 
availability  of  a  suitable  interatomic  potential  is 
crucial  to  a  successful  computer  simulation  study. 

Another  point  concerning  potentials  which  applies  to 


this  work  is  the  difference  between  ideal  and  real  metal 


simulations.  The  properties  of  an  ideal  metal  are  those 

which  pertain  to  metals  in  general.  Real  metal 

properties,  however,  are  those  which  distinguish  a 

particular  metal  from  others.  An  example  of  an  ideal 

metal  property  is  that  divacancies  are  believed  to  be  more 

4 

mobile  than  single  vacancies  in  fee  metals  .  A  molecular 
statics  simulation  which  uses  an  interatomic  potential  for 
a  metal  would  be  expected  to  find  that  the  divacancy 
migration  energy  is  less  than  the  single  vacancy  migration 
energy,  and  hence  that  divacancies  are  more  mobile  than 
single  vacancies.  On  the  the  other  hand,  the  actual 
values  of  the  single  and  divacancy  migration  energies  of  a 
metal  are  real  metal  properties  since  they  are  peculiar  to 
that  metal.  An  interatomic  potential  for  that  particular 
metal  must  be  used  if  one  has  any  hope  of  calculating 
these  values.  This  distinction  between  ideal  and  real 
metal  simulations  is  important  to  this  work  since  we  are 
interested  in  calculating  the  properties  of  real  copper. 

In  this  Chapter  the  potential  functions  cited  and/or 
used  in  this  work  will  be  discussed.  In  the  following 
section  a  general  background  will  be  given  for  empirical 
potentials,  including  some  of  their  strengths  and 
weaknesses.  The  empirical  potentials  are  used  mostly  for 
comparison  here  (some  thermodynamic  calculations  are 
performed)  and  this  discussion  does  not  really  do  them 
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justice.  It  should  also  be  pointed  out  that  "empirical" 

in  this  work  refers  to  potentials  for  which  experimental 

vs.  first  principles  input  is  emphasized,  which  does  not 

31 

follow  standard  usage  ,  and  this  point  will  be  briefly 
discussed.  Following  the  discussion  of  empirical 

potentials,  the  pseudopotential  method  will  be  discussed 
in  a  general  way  and  the  potential  of  Dagens  used  in  this 
work  will  be  described  in  detail.  Some  important 
approximations  that  had  to  be  made  to  use  Dagens' 
potential  are  pointed  out. 

2.2  Empirical  Potentials 

In  this  work,  interatomic  potentials  which  are 

derived  by  a  fit  to  experimental  data  will  be  referred  to 

as  empirical,  although  often  a  distinction  is  made  between 

31 

empirical  and  semi-empirical  potentials  .  In  both  cases, 

the  free  parameters  of  a  functional  form  are  fitted  to 

reproduce  an  experimentally  measured  property  or 

properties.  The  functional  form  of  a  semi-empirical 

potential  is  given  in  advance  of  the  fit  and  is  suggested 

more  by  the  general  properties  of  the  material  (e.g.,  a 

stable  cubic  crystal)  than  by  the  properties  which  are 

actually  being  fitted.  The  functional  form  of  a  true 

empirical  potential,  however,  is  suggested  by  the  property 

31 

or  properties  to  which  it  is  being  fitted  .  The 
potential  consists  of  splines  or  connected  polynomials 
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which  can  represent  a  very  general  function.  Except  for 

32 

the  copper  potential  of  Baskes  and  Melius  briefly 
discussed  below,  the  potentials  discussed  in  this  section 
are  properly  called  semi-empirical  potentials.  Throughout 
the  work,  however,  all  potentials  except  for  the 
interatomic  potential  derived  by  Dagens'  from 
pseudopotential  theory  will  be  called  empirical 
potentials. 

The  general  idea  behind  empirical  potentials  is  to 
assume  that  the  interatomic  potential  can  be  represented 
using  a  fairly  simple  functional  form,  at  least  over 
interatomic  interaction  radii  which  are  important  to  the 
property  or  properties  of  interest.  The  functional  form 
contains  a  few  independent  parameters  which  can  be  fitted 
to  give  the  correct  results  for  a  few  given  properties 
under  very  specific  conditions.  For  example,  one  could 
fit  a  potential  to  the  elastic  constants  of  a  material  at 
a  temperature  of  absolute  zero.  It  is  then  assumed  that 
the  general  form  of  the  potential  is  accurate  enough  to 
use  it  for  other  properties  (interstitial  formation 
energy,  vacancy  migration  energy,  etc.)  and  other 
conditions  (e.g.,  T  >  0). 

One  difficulty  with  empirical  potentials  is  that  the 
functional  parameters  are  not  in  general  unique  and  depend 
on  the  material  property  and  conditions  used  for  the  fit. 
This  problem  can  be  rationalized  in  three  ways.  The  first 


is  that  the  form  of  the  potential  is  not  general  enough  to 
simulate  the  "true"  potential,  i.e.,  the  way  that  real 
atoms  interact.  For  real  materials  this  is  true:  one 
cannot  find  a  single  potential  to  represent  atom-atom 
interactions  under  all  conceivable  conditions.  One 
solution  to  this  could  be  to  generalize  the  potential 
further,  but  simultaneously  fitting  many  properties  is 
difficult  and  it  is  not  clear  that  the  potential  will  be 
more  realistic.  The  second  explanation  is  that  any  form 
of  central  potential  is  insufficient  to  represent 
interatomic  interactions  since  nonlocal  interactions  may 
be  important.  This  means  that  any  form  of  central 
potential,  even  as  general  as  a  table  of  values,  will  not 
be  sufficient  under  all  conditions.  This  is  a 
justification  for  using  empirical  vs.  pseudopotential 
methods.  If  a  single  central  potential  valid  for  all 
properties  does  not  exist,  perhaps  it  is  better  use  an 
empirical  potential  suggested  by  and  fitted  to  the 
properties  in  which  we  are  interested  than  to  go  to  all 
the  trouble  to  derive  an  interatomic  potential  from 
pseudopotential  theory.  Finally,  all  interatomic 
potentials  are  in  the  end  justified  by  their  ability  to 
reproduce  experimental  results.  Even  if  a  potential  has 
some  justification  of  a  first  principles  derivation,  to  be 
useful  it  must  be  able  to  perform  well  in  simulations. 
Hence,  there  are  valid  reasons  for  using  empirical 
potential  in  spite  of  their  disadvantages. 


Among  the  different  forms  of  empirical  potentials. 


three  are  referenced  in  this  work.  The  first  is  the 
Born-Mayer  potential, 


u(r  )  =  L  exp(-yr  )  (2.1) 

•  •  •  • 

1]  ID 

where  r  is  the  interatomic  distance  and  u  is  the  potential 

energy.  This  potential  was  not  used  for  any  calculations 

here,  but  is  cited  for  stacking  fault  results  by  Cotterill 
12 

and  Doyama  .  The  values  of  the  parameters  are 
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L  =  0.053  ev  and  =  13.9  A  ,  from  Gibson  et  ^  .  Since 

the  Born-Mayer  potential  does  not  contain  an  attractive 
part,  an  external  pressure  must  be  applied  to  the  system 
to  maintain  a  given  density.  This  is  a  problem  shared 
with  pseudopotentials,  except  that  the  external  force  can 
in  principle  be  calculated  for  the  pseudopotentials. 

The  second  empirical  potential  cited  in  this  work  in 
the  Morse  potential. 


u(r  )  =  D  {exp[-2<^^r  -r  )]  - 
ij  ij  0 

2  exp[-ct(r  -r  )]}  (2.2) 

ij  0 

where  D,  and  r  are  constants  to  be  specified. 

0 

Cotterill  and  Doyama  used  a  Morse  potential  in  a  study  of 

12 

edge  dislocations  and  stacking  fault  energies  in  copper 
The  potential  was  fitted  to  the  sublimation  energy  of 
copper,  and  the  parameters  for  a  truncation  radius  at  176 


atoms  are  0 


0.3236  ev, 


-1 

=  1.2941  A  ,  and 


r  *  2.9133  A. 

0 

The  third  empirical  potential  is  a  modified  Morse 

29 

potential  used  by  Schober  ^  al  to  calculate 
interstitial  properties.  The  functional  form  of  this 
potential  is 


u(r  )  =  D  {exp[-2oi<r  -r  )]  - 
ij  ij  0 

2  exp[-dc(r  -r  )]}  + 
ij  0 

4 

B(r  -r  )  H(r  -r  )  +  V  (2.3) 

ij  NN  NN  ij  0 

where 


K(x)  = 


X  <  0 

X  >  0 


(2.4) 


The  extra  terms,  compared  with  equation  (2.2),  makes  the 

potential  softer  in  the  nearest  neighbor  region  and 

provides  two  more  parameters  which  can  be  used  for  fitting 
29 

properties  .  This  potential  was  fitted  to  the  vacancy 

formation  energy  and  elastic  constants  of  copper  and  the 

potential  parameters  are  D  =  0.227  ev,  r  =  0.71555a, 

0 


QL  = 


r 

NN 

this 
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7.6499a  ,  B  =  -1200ev  /  a  ,  V  =  0.0233  ev,  and 

0 

1/2 

a  /  2  ,  where  a  is  the  lattice  constant.  Since 

potential  was  fitted  to  more  recent  values  of  the 
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vacancy  formation  energy  and  was  also  fitted  to  the  bulk 

modulus,  it  is  used  in  Chapter  5  to  do  thermodynamic 

calculations  of  copper  for  comparison  with  the  Dagens' 

interatomic  potential  discussed  below. 

An  empirical  potential  recently  derived  by  Baskes  and 
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Melius  was  also  considered  .  Baskes  and  Melius  did  a 

true  empirical  fit  for  seven  fee  metals  to  a  large  number 

of  properties,  and  had  a  separate  term  for  the  cohesive 

electron  energy  (the  structure- independent  term  in  the 

pseudopotential).  Unfortunately,  they  apparently  used 

outdated  values  for  the  vacancy  formation  and  migration 

4 

energies  in  the  fit  for  copper  (see  Balluffi  and 

5 

Wienhold  for  recent  estimates  of  these  values).  De  Leeuw 

et  al  used  this  potential  to  calculate  copper  properties 

near  the  melting  point  and  found  that  a  significant 

external  pressure  had  to  be  applied  to  maintain  the 

34 

correct  density  .  Because  of  these  problems  this 

potential  was  not  used  in  this  work. 


2.3  Pseudopotentials 


There  is  a  long  history  behind  the  pseudopotential 

method,  and  the  classic  work  was  published  by  Harrison  in 
16 

1966  .  The  pseudopotential  method  is  an  attempt  to 


derive  interatomic  potentials  from  first  principles.  It 
takes  advantage  of  the  f ree-electron  nature  of  the 
conduction  electrons  in  metals,  which  are  free  to  roam 


throughout  the  metal  instead  of  being  associated  with 

individual  atoms.  The  metal  is  separated  into  two 

components:  the  conduction  electrons,  which  are 

represented  by  some  type  of  average  interaction,  and  the 

16 

metallic  ions,  whose  electrons  are  tightly  bound  .  The 

strong  potential  inside  the  metallic  ions  is  then  replaced 

by  a  weaker  potential  which  preserves  the  energy 
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eigenvalues  of  the  conduction  electrons  .  The  original 

pseudopotential  method  proposed  by  Harrison  did  not 

require  experimental  information  to  be  introduced  into  the 

derivation  of  potential.  Since  the  late  1960's,  however, 

modifications  of  the  method  which  incorporate  some 

experimental  data  have  proven  more  successful  for  some 

35,  36,  37 

metals  than  the  pure  pseudopotential  approach 

After  discussing  some  general  aspects  of  pseudopotentials, 

the  potential  used  in  this  work  will  be  covered. 

Interatomic  potentials  derived  from  pseudopotential 

methods  have  some  common  characteristics.  The  first  is 

that  the  potentials  are  oscillatory  beyond  a  range  of 

38 

about  one  lattice  constant  .  These  oscillations  converge 

to  the  Friedel  oscillations  at  long  range.  Whether  or  not 

these  oscillations  are  physically  significant  is  not 

39,  37 

well-established  .  One  practical  consequence  of  the 

oscillations  is  that  these  potentials  are  long-ranged  and 
must  be  truncated  to  be  used  for  molecular  dynamics  and 
statics.  This  can  lead  to  convergence  difficulties. 


especially  for  properties  which  depend  on  the  long-range 
part  of  the  potential. 

A  second  characteristic  of  the  potentials  is  they 
contain  structure-dependent  and  independent  components 
which  arise  from  the  treatment  of  the  free  electrons.  If 
the  total  cohesive  energy  of  the  crystal  is  included  in 
the  interatomic  pair  potential  u  ,  the  total  energy  E  can 

ij 

be  written 

E  =  23  u(r  )  (2.5) 

i<j  ij 

where  r  is  the  distance  between  atoms  i  and  j.  Equation 

ij 

(2.5)  is  valid  for  many  empirical  pair  potentials 

including  the  Morse  and  modi f ied-Morse  potentials 

discussed  in  the  previous  section.  In  metals,  however, 

and  particularly  when  the  pseudopotential  method  is  used, 

40 

the  total  cohesive  energy  is  given  by 

E  =  u(r  ;V)  +  U(V)  (2.6) 

i<j  ij 

The  first  term  in  (2.6),  similar  to  the  right-hand  side  of 

(2.5) ,  is  the  structure-dependent  component  of  the  total 

energy  and  depends  on  the  positions  of  the  atoms. 
This  part  of  the  potential  is  used  directly  in  molecular 
dynamics  since  it  describes  how  one  atom  interacts  with 
another.  The  second  term  of  (2.6)  is  the  structure- 

independent  component  and  depends  only  on  the  system 
volume  and  not  on  the  positions  of  the  atoms.  It 


physically  represents  the  cohesive  energy  provided  by  the 

conduction  electrons,  which  are  not  directly  simulated. 

It  is  not  necessary  to  know  the  structure-independent  part 

of  the  potential  to  calculate  many  properties  of  the 

material.  For  example,  if  we  wanted  to  determine  the 

stable  lattice  configuration  or  calculate  defect  formation 

energies  at  a  given  density,  we  would  only  need  the 

structure-dependent  part  of  the  potential.  For  this 

reason  the  structure- independent  component  is  not  always 

provided  in  derivations  of  pseudopotentials.  This  is 

unfortunate  since  the  volume  derivative  of  the  total 

energy  is  the  system  pressure,  and  this  derivative 

obviously  depends  on  the  second  term  of  (2.6).  When  the 

structure-independent  component  is  known,  it  adds  an  extra 

complication  to  direct  thermodynamic  and  mechanical 

property  calculations.  When  it  is  unknown,  some 

thermodynamic  properties  simply  cannot  be  evaluated. 

To  make  matters  worse,  the  structure-dependent 

component  of  the  potential  is  also  volume  dependent. 

Hence,  the  form  of  the  interatomic  potential  u(r;  V) 

depends  on  the  density  of  the  system.  If  simulations  are 

done  at  the  same  volume  V  ,  a  single  interatomic  potential 

0 

u(r;  V  )  can  be  used  and  the  the  structure-independent 
0 

component  U(V  )  is  not  necessary  as  long  as  it  is 
0 

sufficient  to  calculate  the  total  energy  within  a 
constant.  To  compare  calculations  at  a  different  volume 
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V  ,  however,  we  not  only  need  to  know  the  structure- 
1 

independent  component  U(V  ),  but  we  also  need  a  different 

1 

interatomic  potential  u(r;  V  ). 

1 

Another  complication  in  using  pseudopotentials  in 
molecular  dynamics  and  molecular  statics  calculations  is 
that  they  generally  do  not  come  in  a  simple  functional 
form.  The  interatomic  potential  calculated  using  the 
pseudopotential  method  consists  of  a  table  of  values  of 
energy  vs.  interatomic  radius.  In  principle  this  is  not  a 
problem;  one  can  calculate  the  table  to  any  desired 
accuracy  with  sufficient  computer  resources,  and  the  table 
can  be  interpolated  using  standard  numerical  methods  such 
as  splines.  However,  it  requires  a  considerable  amount  of 
space  to  publish  a  table  of  sufficient  accuracy  for  all 
calculations,  and  programming  would  be  tedious  and  prone 
to  error  using  such  a  table  from  a  paper.  The 

pseudopotential  method  is  difficult  enough  to  preclude 
each  worker  from  repeating  the  calculations  to  get  the 
interatomic  potential  function,  especially  using  the 
advanced  methods  necessary  to  get  the  best  interatomic 
potential.  One  solution  to  this  dilemma  is  to  fit  the 
interatomic  potential  generated  by  the  pseudopotential 
method  to  a  very  general  functional  form  suggested  by  the 
physics  of  interatomic  interactions  and  the  structure  of 
the  interatomic  potential  itself.  Of  course,  the  function 
to  which  we  are  fitting  the  pseudopotential  interatomic 
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potential  should  not  be  so  complicated  that  it  also  causes 
difficulties. 


In  spite  of  all  of  these  problems  there  still  is  much 
interest  in  developing  and  using  interatomic  potentials 
derived  from  pseudopotential  theory.  The  reasons  are  that 
the  potentials  are  better  justified  theoretically  than 
empirical  potentials,  and  also  that  the  pseudopotential 
methods  provide  a  theory  of  interatomic  interactions  in 
metals. 

The  interatomic  potential  used  in  this  work  is  based 

13,  14,  15 

on  the  work  of  Oagens  .  This  was  the  culmination 

of  attempts  to  account  for  the  d-band  in  the  noble 
41,  42 

metals  .  Harrison  had  previously  noted  that  a  basic 

assumption  behind  the  pseudopotential  method  was  that  the 

core  and  conduction  electrons  could  be  cleanly  separated 

such  that  the  core  of  the  atom  could  be  treated  as  an 
43 

ion  .  The  assumption  only  strictly  holds  for  the  alkali 

and  polyvalent  metals  for  which  reliable  interatomic 

potentials  have  been  derived  using  the  pseudopotential 

method.  Unfortunately,  the  assumption  breaks  down  in 

noble  and  transition  metals. 

The  results  of  Dagens'  calculations  were  fitted  to  a 

1 

general  interatomic  potential  of  the  form 


u(r)  =  {1  -  exp[-/X(r-r  )  ]}  x 

1 


C  \  S  sin^ 

1  \  cos^  1 

C  +  —  -  +  -  + 

0  2  3  5 


B  exp(-jgr  ) 

- 

n  3-n 

r 

D  exp(-iir) 
Aer«  76  =  Y  +  2k^r. 


(2.7) 


The  parameters  for  the  fit  are  given  in  Table  2-1.  The 

parameters  for  silver  are  also  included  in  Table  2-1,  as 

this  potential  will  be  used  in  Chapter  8.  The  copper 

potentials  discussed  in  this  Chapter  are  plotted  in  Figure 

2-1.  The  differences  between  the  various  potentials  is 

striking  when  it  it  considered  that  they  are  all  for  the 

1 

same  metal.  As  Lam  points  out,  the  copper  potential  has 

a  minimum  near  the  first  nearest  neighbor  radius.  Whether 

or  not  this  is  physically  realistic  is  not 

37 

well-established 

It  was  not  possible  to  determine  the  structure- 

independent  part  of  the  potential  for  this  work; 

44 

apparently,  this  is  a  very  difficult  task  .  This  does 
not  affect  the  defect  calculations  presented  in  this  work. 
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16.55 


17.669 


B 

0 


B 

-0.0033 

-0.459 

1 

B 

0.0825 

-0.2044 

2 

B 

-0.1576 

-0.02728 

3 

0.74 

oo 

o 

D 

0.1246 

0.2819 

k 

0.7189 

0.6387 

F 

Table  2-1:  Interatomic  Potential  Parameters 

1 

for  Noble  Metals  (from  Lam  et  al  ) 


a 


but  it  made  it  impossible  to  calculate  the  pressure  for 
the  thermodynamic  calculations.  This  problem  will  be 
discussed  in  more  detail  in  Chapter  5. 

As  pointed  out  earlier  in  this  Chapter,  the 

interatomic  potentials  derived  using  the  pseudopotential 

method  are  long-ranged.  To  overcome  this  difficulty, 

45  .  r 

Duesbury  et  al  developed  a  damping  factor  for  the 

potential  u(r)  such  that 


eff 

u  (r)  =  u 

2  2 

(r)  exp(-d  r  ) 

(2.8) 

eff 

where  u  (r)  is 

the  potential 

actually  used  in 

1 

the 

molecular  statics 

calculations. 

Following  Lam  , 

the 

copper  calculations 

using  the  Dagens'  potential  were 

done 

using  d  =  0.25  a  ,  where  a  is  the  lattice  constant. 


Chapter  2 

Molecular  Dynamics 


3.1  Introduction 

The  thermodynamic  properties  of  copper  were 

calculated  in  this  work  using  a  technique  called 

computational  molecular  dynamics  (MD) .  In  MD  one  solves 

the  equations  of  motion  of  N  atoms  subject  to  a  given 

interatomic  potential  using  a  computer.  Given  N  atoms 

with  coordinates  r  and  velocities  v  ,  1  5  i  £  Nf 

i  i 

solves  the  set  of  equations 


dv  (t) 

^  T 

m  -  =  -  Z- 


du(r  ) 

ij 


(3.1) 


dt 


j?^i  dr 


ID 


for  r  (t)  and  v  (t),  where  m  is  the  atom  mass,  r  is  the 
i  i  ij 

distance  between  atoms  i  and  j,  and  u(r)  is  the 


interatomic  potential.  Many  properties  of  the  system  such 
as  its  temperature,  internal  energy,  and  pressure,  can  be 


calculated  since  they  can  be  expressed  in  terms  of  r  (t) 

i 

and  V  ( t ) . 
i 

The  system  that  one  simulates  using  conventional  MD 


has  a  constant  number  of  atoms  N,  a  constant  volume  V,  and 
a  constant  total  energy  E.  This  type  of  system  is  called  a 


microcanonical  ensemble.  The  ensemble  itself,  which 
consists  of  all  the  possible  configurations  of  the  system 
with  the  given  N,  V,  and  E,  is  uniquely  specified  by  these 
values.  Each  configuration  of  the  system  has  it's  own 
temperature  and  pressure,  however,  and  the  temperature  and 
pressure  of  the  ensemble  is  the  average  of  these  values 
over  all  the  possible  configurations  in  the  ensemble. 
Even  though  the  temperature  and  pressure  of  the  ensemble 
are  well-defined  and  unique,  one  cannot  reliably  estimate 
their  values  without  taking  an  average  over  many 
configurations  or  MD  time  steps.  Consequently,  one  cannot 
know  the  temperature  and  pressure  of  the  ensemble  which  is 
being  simulated  until  after  the  run  is  completed. 

These  characteristics  of  conventional  MD  affect  its 
use  for  thermodynamic  property  calculations.  As  pointed 
out  above,  it  is  necessary  to  average  the  temperature  and 
pressure  (and  other  properties  which  fluctuate)  over  many 
time  steps  to  determine  the  ensemble  averages.  There  is 
no  practical  way  to  average  over  all  possible 
configurations,  since  (classically)  there  are  an  infinite 
number  of  them,  and  the  best  one  can  do  is  get  an  estimate 
of  a  fluctuating  property.  It  is  important  that 
"unlikely"  configurations,  i.e.,  configurations  which  are 
very  unlikely  to  occur,  are  excluded  because  their 
property  values  will  distort  the  averages.  These  unlikely 
configurations  have  the  greatest  chance  of  occurring  at 


the  very  beginning  of  the  simulation,  since  as  the 
simulation  proceeds  it  will  approach  more  probable 
configurations.  Hence,  it  is  important  to  select  the 
initial  configuration  properly.  The  consequence  of  not 
knowing  the  temperature  and  pressure  of  the  ensemble  one 
is  simulating  until  well  into  the  simulation  is  that  is 
difficult  calculate  ensembles  with  a  given  pressure  and 
temperature.  This  is  unfortunate  since  experimental 
properties  of  a  material,  especially  those  for  solids,  are 
usually  presented  at  a  given  temperature  and  for 
atmospheric  pressure.  Without  special  techniques,  the 
best  one  can  expect  to  do  is  to  calculate  ensembles  for 
various  values  of  V  and  E  and  interpolate  to  the  desired 
temperature  and  pressure.  In  this  work,  velocity-scaling 
is  used  to  set  the  simulation  temperature  and  the  flexible 
boundary  technique  is  used  to  set  the  simulation  pressure. 

In  this  Chapter  the  molecular  dynamics  method  will  be 

discussed.  The  integration  scheme  used  to  solve  (3.1)  for 

r  (t)  and  v  (t)  is  presented  first.  This  is  followed  by  a 
i  i 

discussion  of  the  initial  and  boundary  conditions,  which 
includes  the  selection  of  the  initial  configuration  and 
velocity-scaling.  The  equations  for  calculating  basic 
properties  are  given,  and  computational  time-saving 
methods  used  in  this  work  are  discussed.  Finally,  the 
flexible  border  method  is  outlined. 


A 


3.2  Integration  Scheme 


The  fundamental  feature  of  molecular  dynamics  is  the 

solution  of  the  equations  of  motion  given  by  (3.1). 

Except  for  simple  cases  it  is  not  possible  to  solve  these 

equations  analytically  and  it  is  necessary  to  resort  to 

numerical  methods.  Since  the  solution  of  these  equations 

is  so  basic  to  computational  molecular  dynamics,  the 

selection  of  the  numerical  integrating  scheme  is 

important.  A  detailed  study  of  integration  methods  for  MD 

46 

applications  has  been  made  by  Beeman  .  The  method  which 

appears  to  be  is  a  Gear-Nordsieck  method  discussed  by 
47 

Haile  .  This  method  has  three  basic  steps: 

1.  Prediction  -  use  values  of  positions  and 
derivatives  at  time  t  to  estimate  the  positions 
and  derivatives  and  t+^t. 

2.  Evaluation  -  evaluate  the  potential  function, 
i.e.,  the  right  hand  side  of  (3.1),  at  the 
predicted  positions. 

3.  Correction  -  using  the  results  of  the  evaluated 
potential,  correct  the  positions  and 
derivatives. 


(n)  th 

Let  r  (t)  be  the  n  derivative  of  the  position  of 

P(n) 

atom  i  and  time  t,  and  let  r  (t+^t)  be  the  predicted 

i 

th 

value  of  the  n  derivative  of  the  position  of  atom  i  at 

time  t+At.  For  an  integrator  of  order  p,  we  need  the 

P(n) 

derivatives  up  to  n  *  p.  We  estimate  r  (t+^t)  by 

i 

expanding  the  Taylor  series  in  r  (t)  about  t,  where 


(t)  (/^t) 


where  p  is  the  order  of  the  integrator  we  wish  to  use. 

The  next  step,  and  the  core  of  a  molecular  dynamics 

P(n) 

code,  is  the  evaluation  of  the  potential  for  r  (t+^t). 

i 

Most  of  the  computer  time  in  a  molecular  dynamics  run  is 

spent  in  this  step,  and  many  of  the  details  of  the 

implementation  are  buried  here.  These  details,  such  as 

spatial  boundaries  which  are  imposed  on  the  system  and 

time-saving  techniques  for  the  evaluation  of  the 

potential,  will  be  discussed  in  subsequent  sections. 

The  evaluation  of  the  right  hand  side  of  (3.1) 

yields,  after  dividing  by  the  mass  m,  the  second 
(2) 

derivatives  r  (t+^t).  The  difference  between  predicted 

P(2)  ^  (2) 

value  r  (t+^t)  and  the  calculated  value  r  (t+^t) 

i  i 

(2) 

yields  an  estimate  of  the  error  in  the  prediction  ^  , 

where 


(2)  (2)  P(2) 

AL  “I  (t+£jt)  -  r  (t+^t) 
i  i 


(3.3) 


The  error  ^  is  now  used  to  make  corrections  to  all  the 


derivatives  where 


r  (t+2^t)  =  r  (t+^t)  + 


Za?: 

n 


n-2 

(2)  (^t)  n! 


(3.4) 


where  cl  is  given  in  Table  3-1.  Using  the  values  r  we 
i  i 

now  make  another  time  step  by  repeating  the  process. 

cc  Order  p 


19/120 


1/12 


3/16 

251/360 


11/18 


1/60 


Table  3-1:  Values  of  ot  for  Gear  Predictor-Corrector 


Algorithm 


We  found  it  useful  in  some  phases  of  this  work  to 
extend  the  Gear  algorithm  to  perform  in  a  variable 
time-stepping  mode,  where  the  integrator  varies  the  size 
of  the  time-step  automatically.  The  idea  behind  variable 
stepsize  integrating  algorithms  is  to  estimate  the  error  E 
in  making  a  timestep  by  comparing  the  results  of  two 


integrators  of  orders  q  and  q+1.  If  the  estimated  error  E 

is  larger  than  the  maximum  allowed  error  E  ,  we  repeat 

max 

the  step  at  a  smaller  stepsize.  If  the  error  is  less  the 

maximum  allowed,  we  accept  the  step  and  change  the 

[q] 

stepsize  based  on  the  difference.  Let  r  (t+^t)  and 
[q+1] 

r  (t+^t)  be  the  positions  calculated  using  integrators 
i 

of  order  q  and  q+1.  The  estimated  error  for  the  time  step 
is  given  by 


E  -  max 
i,k 


[q+1]  [q] 

r  (t+^t)  -  r  (t+^t) 
i  ,k  i  ,k 

th 


(3.5) 


where  r  is  the  k  component  of  r  .  The  new  stepsize 
i  ,k  i 
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^t  is  calculated  from  the  old  stepsize  ^t  by 
new  old 


0.1  ^t 


old 


At 


^  At 


new 


old 

2.0  ^t 


old 


H  1  0.1 
0.1  <  <  2.0 
2.0  <  X 


(3.6) 


and 


l/(q+l) 

=  0.9  (E  /  E) 
max 

The  new  stepsize  ^t  is  used  if  the  maximum  error  is 

new 

exceeded  and  we  have  to  repeat  the  step,  or  if  we  don't 


exceed  the  maximum  and  the  step  is  accepted. 


The  usefulness  of  the  variable  stepsize  method  is 
limited  in  the  applications  in  this  work.  The  major 
drawback  of  the  method  is  that  two  function  evaluations 
are  necessary,  which  automatically  makes  it  almost  twice 
as  slow  as  the  fixed  stepsize  method.  The  time  step  size 
using  the  variable  step  algorithm  is  fairly  constant 
during  thermodynamic  property  calculations,  and  certainly 
does  not  vary  by  a  factor  of  two.  The  advantage  of  using 
a  variable  stepsize  method  is  that  one  does  not  have  to 
select  a  stepsize  and  hence  the  amount  of  work  the  user 
has  to  do  and  the  chance  of  user  error  are  decreased.  We 
took  advantage  of  this  in  some  of  the  property 
calculations;  the  variable  stepsize  integrator  was  used 
for  25-100  steps,  and  the  smallest  stepsize  necessary 
during  that  period  was  used  for  the  remainder  of  the  run 
in  a  fixed  stepsize  mode. 


3.3  Boundary  and  Initial  Conditions 


The  specification  of  a  differential  equation  problem 
is  incomplete  without  boundary  and  initial  conditions. 
Equation  (3.1)  is  a  second  order  ordinary  differential 
equation  in  time  and  requires  2Nd  initial  conditions, 
where  N  is  the  number  of  atoms  and  d  is  the  number  of 
dimensions.  Physically,  these  initial  conditions  are  the 
initial  positions  and  velocities  of  the  atoms. 

As  pointed  out  in  the  Introduction  to  this  Chapter, 


it  is  essential  to  exclude  "unlikely"  configurations  when 

taking  property  averages.  A  "likely"  configuration  is  one 

in  which  the  atoms  are  evenly  distributed  in  the  system 

volume  and  the  velocities  are  distributed  in  a  Boltzmann 

distribution.  It  is  possible  to  set  up  a  configuration 

such  as  this  by  assigning  the  atoms  to  lattice  sites  of  a 

lattice  type  such  as  fee,  and  assigning  initial  velocities 

randomly  according  to  a  Boltzmann  distribution.  At  a 

nonzero  temperature,  however,  one  does  not  expect  to  find 

the  atoms  located  at  perfect  fee  sites  but  slightly 

displaced  from  these  positions.  The  practical  solution  to 

this  problem  is  to  note  that  as  the  simulation  proceeds, 

one  expects  that  the  system  will  approach  more  likely 

configurations,  and  hence  regardless  of  the  initial 

configuration  the  system  eventually  approaches 

configurations  which  are  suitable  for  property  averages. 

This  process  is  called  "equilibration"  and  may  take  a  few 

47 

hundred  time  steps  .  For  some  property  calculations  the 
atoms  were  set  at  fee  lattice  sites  and  the  velocities 
were  assigned  from  a  uniform  distribution.  The  simulation 
was  started,  but  property  averages  were  not  taken  until  a 
few  hundred  time  steps  had  passed.  In  other  calculations, 
the  final  configuration  from  previous  simulations  was  used 
as  the  initial  configuration. 

So  far  no  mention  has  been  made  about  the  assignment 
of  the  initial  velocities  of  the  system  except  to  say  that 


they  are  selected  from  some  type  of  distribution  or  from 

previous  runs.  An  important  characteristic  of  MD  using 

the  periodic  boundaries  discussed  below  is  that  linear 

momentum  is  conserved.  The  initial  linear  momentum 

introduced  into  the  system  is  maintained  throughout  the 

simulation.  If  this  momentum  is  nonzero,  the  total  system 

will  translate  over  time  and  the  temperature  will  be 

calculated  incorrectly  since  it  will  include  the 

translational  velocity.  Hence,  it  is  important  to  insure 

that  the  linear  momentum  of  the  system  is  zero. 

The  initial  velocities  of  the  system  provide  a  way  of 

setting  the  approximate  temperature  of  the  system. 

Although  it  is  not  possible  before  actually  doing  the 

calculations  to  know  the  final  average  temperature  of  an 

ensemble,  it  is  possible  to  calculate  and/or  set  the 

instantaneous  temperature  of  a  configuration.  If  T  is  the 

current  configuration  temperature  and  T  is  the  desired 

des 

temperature,  one  can  set  the  instantaneous  temperature  to 

T  by  multiplying  all  of  the  atom  velocities  by 
des 
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(T  /T)  .  This  is  called  velocity-scaling  .  In  this 

des 

work,  velocity-scaling  was  used  over  a  few  hundred 
preliminary  time  steps,  during  which  the  system 
equilibrated  and  no  property  averages  were  calculated. 
The  system  was  then  allowed  to  rest  for  twenty  or  more 
steps,  and  then  property  calculations  were  started.  The 


final  average  temperature  was  usually  within  ten  percent 
of  the  desired  temperature. 


Since  there  are  no  spatial  derivatives  in  (3.1) 

(recall  that  u(r  )  is  a  known  function),  it  would  seem 

i 

that  spatial  boundaries  should  not  be  a  problem.  The 

problem  arises,  however,  because  we  are  trying  to  simulate 

a  bulk  material  containing  a  large  number  of  atoms,  but 

are  computationally  limited  to  a  few  hundred  atoms  and 

simulation  volumes  with  edge  lengths  of  a  few  lattice 

constants.  On  a  scale  this  size,  where  almost  every  atom 

can  interact  across  the  simulation  boundary,  the  boundary 

conditions  are  very  important.  To  overcome  this 

difficulty,  periodic  boundaries  in  which  the  simulation 

volume  is  surrounded  by  images  of  itself  as  shown  in 

47 

Figure  3-1  are  commonly  employed  .  In  two  dimensions, 
every  atom  has  eight  images.  An  atom  interacts  with  the 
nearest  image  of  another  atom.  This  is  shown  in  Figure 
3-1  where  atom  A  interacts  with  the  image  of  B  in  the 
lower  left  volume. 

3.4  Time-Saving  Methods 

It  is  possible  to  save  considerable  amounts  of 
computer  time  by  taking  advantage  of  some  time  saving 
techniques.  In  this  section  we  will  discuss  the 

techniques  used  in  this  code. 
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Figure  3-1:  Periodic  Boundaries 


The  simulation  cell  in  the  center  is  surrounded 
by  images  of  itself.  Atom  A,  only  shown  in  the 
center  cell,  does  not  interact  directly  with  atom 
B  but  with  one  of  the  images  of  B. 


3.4.1  Verlet  Neighbor  List 

As  shown  in  Equation  (3.1),  it  is  necessary  to 
calculate  the  force  between  each  atom  during  every  time 
step.  In  principle  this  requires  N(N-l)/2  evaluations  of 
the  potential  function.  As  noted  in  Chapter  2,  however, 
the  potential  function  is  normally  truncated  so  that  there 
are  atoms  in  the  simulation  which  do  not  directly  interact 
with  one  another  during  some  part  of  the  run.  The 
distance  between  the  closest  images  of  two  atoms  was 
greater  than  the  cutoff  distance  for  the  potential,  it  is 
evident  that  the  potential  between  the  two  atoms  is 
trivial  calculate  since  it  is  zero.  At  this  point, 

however,  computer  time  has  already  been  wasted  calculating 
the  distance  between  two  noninteracting  atoms.  Since  we 
have  N(N-])/2  interactions  to  calculate,  this  wastage  can 
be  considerable. 

Of  the  suggested  solutions  to  this  problem,  the 
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standard  Verlet  list  was  used  in  this  work  .  The  Verlet 
neighbor  list  is  used  to  keep  track  of  the  atoms  which 

have  a  nonzero  interaction.  One  can  think  of  the  Verlet 

list  as  a  NxN  truth  matrix  T  ,  where  T  =  1  or  0  if 

ij  ij 

atoms  i  and  j  do  or  do  not  interact,  respectively.  T  is 

ij 

evaluated  by  actually  calculating  and  checking  the 
distances  between  all  the  atoms.  For  a  subsequent  given 
number  of  time  steps,  instead  of  calculating  the  distance 


IS  1  or 


between  atoms  i  and  j,  we  just  check  to  see  if  T 

ij 

0,  which  is  considerably  quicker  (and  only  calculate  the 
distance  if  it  is  necessary  to  evaluate  the  potential). 
Since  we  are  doing  a  dynamical  calculation,  it  is  evident 
that  T  will  eventually  become  out-of-date  (in  the  worst 

ij 

case,  after  one  time  step).  We  overcome  this  difficulty 

by  including  atoms  from  slightly  larger  radii  than  the 

potential  cutoff  radius.  Although  we  have  a  few  extra 

atoms  in  the  list,  we  still  calculate  fewer  than  N(N-l)/2 

radii  and  the  list  will  remain  up  to  date  for  a  greater 

number  of  time  steps.  In  this  work,  the  cutoff  range  for 

the  Verlet  list  was  set  at  l.lr  ,  where  r  is  the  cutoff 

C  C 

radius  for  the  potential,  and  the  list  was  updated  every 
ten  potential  evaluations.  It  is  probably  unnecessary  to 
to  update  the  list  this  often  since  we  were  usually 
simulating  solids,  but  it  was  found  that  the  percentage  of 
time  spent  updating  the  list  at  this  frequency  was 
relatively  small. 

The  truth  matrix  T  is  not  usually  implemented  as  a 

ij 

matrix  on  the  computer  for  a  number  of  reasons.  Since  T 

ij 

is  symmetric  and  the  diagonal  is  all  zeros,  we  actually 
only  need  to  know  the  strictly  upper  triangle  of  T  .  If 

ij 

the  potential  cutoff  radius  is  much  smaller  than  the 
simulation  volume  size,  many  elements  T  are  0  and  the 

ij 

matrix  is  sparse.  Even  though  it  is  much  quicker  to  test 
T  than  calculate  and  test  the  distance  be  atoms  i  and  j. 


we  still  have  to  test  (N-l)(N-2)/2  values.  Finally,  the 

2 

niimber  of  elements  in  T  is  proportional  to  N  ,  so  the 

ij 

matrix  quickly  gets  large  as  N  increases. 

The  standard  way  of  implementing  the  neighbor  list  is 

illustrated  schematically  in  Figure  3-2.  The  truth  table 

T  is  shown  for  a  hypothetical  system  of  N  =  5  atoms.  As 
i  j 

pointed  out  above,  the  upper  triangle  of  T  is  all  we 

ij 

need  to  know.  The  idea  behind  the  method  is  to  keep  track 
of  only  the  I's  in  T  .  This  is  done  with  five  vectors, 

ij 

one  for  each  atom,  where  we  note  that  the  last  vector  is 
always  null.  The  final  step  is  to  pack  the  vectors  in  a 
sequential  list,  where  we  keep  track  of  the  start  of  an 
atom's  list  using  another  vector. 

3.4,2  Evaluation  of  the  Potential 

1 

The  potential  function  given  by  given  by  Lam  and 
discussed  in  Chapter  2  contain  many  transcendental 
functions  which  are  time-consuming  to  calculate  on  a 
computer.  Since  most  of  the  computer  time  in  MD  is 
usually  spent  calculating  the  potential  for  each  atom-atom 
interaction  for  each  time  step,  the  potential  in  its 
analytic  form  is  a  serious  practical  problem.  The  common 
solution  is  to  take  the  complicated  potential  and  fit  it 
to  a  functional  form  which  is  computationally  easy  to 
evaluate,  such  as  polynomials,  splines,  or  table-lookup. 
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We  used  a  spline  fit  which  had  maximum  relative  error  of 
-4 

less  than  10  over  the  radii  of  interest. 

In  addition  to  using  the  spline  fit  to  save  time,  ve 
also  used  another  computer-related  method.  During  each 
time  step  in  an  MD  simulation  it  is  necessary  to  calculate 
to  force  between  each  atom.  The  general  procedure  is  to 
calculate  the  radius  r  between  the  two  atoms  and  then 
calculate  the  force  for  radius  r.  This  requires 
calculating  this  each  component  of  r  and  then  finding  the 
Euclidean  norm  by  taking  the  square  root  of  the  sum  of  the 
squares  of  the  components.  Depending  on  the  computer 
used,  the  step  that  often  takes  the  longest  is  taking  the 
square  root.  If  the  potential  function  is  particularly 
simple  the  calculation  of  r  may  take  longer  than 
calculating  the  force.  For  this  reason  it  is  possible  to 
save  a  considerable  amount  of  computer  time  if  the 
evaluation  of  the  square  root  can  be  avoided. 

Let  U(r)  be  a  potential  function,  where  r  is  the 
radius  between  two  particles.  During  an  MD  run  we  need  to 
calculate 


ij  dU 

F  (r  )  -  -  —  ij  X 
k  dr  r=r 

i  j  th 

where  r  is  the  k  component 
k 

another  potential  W(s) 


ij 

r 

k 

ij 

r 

i  j 

of  r  -r  . 


(3.7) 

Also  define 


t 


: 


V' 
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W(s)  =  U(r)  (3.8) 


2 

where  s  =  r  .  It  is  readily  shown 


dW(s) 

1 

dU{r) 

(3.9) 

—  A 

ds 

2r 

dr 

and 


i j  dW ( s )  i j 

F  (r  )  =  -2  -  r 

(3.10) 

k  ds  k 

In  most  simulations  involving  thermodynamic  property 
calculations  the  virial  term  must  be  evaluated  to 
calculate  the  pressure.  This  can  be  expressed 

dU(r)  dW(s) 

r  -  =  2s  - 

dr  ds 

Hence,  if  the  potential  function  can  be  expressed  as  a 

2 

function  of  s  =  r  ,  we  can  calculate  the  energy,  force 

components,  and  virial  term  without  taking  a  square  root. 

In  this  work  this  involved  fitting  the  potential  to  s 

instead  of  r.  The  run  times  on  the  Compupro  8086  system 

were  reduced  by  approximately  forty  percent.  We  also  note 

that  the  commonly  used  Lennard-Jones  6-12  potential  is 

2 

easily  expressed  as  a  function  of  r  ,  although  the  shifted 

force  method  in  its  common  form  cannot  since  that  puts  an 

47 

r  term  in  the  potential  .  It  is  possible  to  shift  the 

2 

force  using  an  r  term,  however,  which  probably  has  as 
much  physical  significance  as  the  standard  way. 


I 
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3.5  Thermodynamic  Property  Calculations 

Three  basic  thermodynamic  properties  are  normally 

calculated  in  MD  simulations:  the  total  energy  E,  the 

temperature  T,  and  the  pressure  P.  In  the  microcanonical 

simulations  performed  using  conventional  molecular 

dynamics,  the  simulation  volxwie  is  fixed  and  hence  the 

mass  density  is  predetermined.  Since  the  total  energy 

should  be  constant,  it  is  also  a  useful  measure  of  the 

correctness  of  the  computer  code.  In  this  section  we  will 

outline  how  E,  T,  are  P  are  calculated  using  molecular 

47 

dynamics  following  Haile  .  Other  properties  which  are 
also  calculated  simultaneously  with  the  thermodynamics 
properties,  the  radial  distribution  function,  order 
parameter,  and  the  mean-square  displacement,  will  also  be 
discussed  in  this  Section. 

The  total  energy  E  is  the  sum  of  the  kinetic  energy 
K.E.  and  potential  energy  P.E.  of  the  system, 

E  =  K.E.  +  P.E.  (3.11) 

The  instantaneous  kinetic  energy  is  given  by 

3m  2 

K.E.{t)  =  — 2.V  (t),  (3.12) 

2  i  i 

i.e.,  the  sum  of  the  kinetic  energies  of  each  atom.  The 


instantaneous  potential  energy  is  the  sum  of  the 
interactions  between  each  atom, 


P.E.(t)  *21  u(r  ) 
j<i  ij 


(3.13) 


Since  we  are  solving  Newton’s  equations  of  motion,  we 

expect  the  total  energy  E  to  be  constant,  but  because  of 

-4 

numerical  errors,  E  varies  by  about  10 

The  instantaneous  temperature  T  can  be  extracted  from 
the  kinetic  energy  since 


3 

-  NkT  =  K.E.  (3.14) 

2 

and  solving  for  T 


2 

T  =  -  K.E.  (3.15) 

3Nk 

The  simulation  pressure  P  is  more  difficult  to 
calculate  than  U  and  T.  The  instantaneous  pressure  P  is 
given  by 


du(r  ) 

NkT 

V  3V  i<j  ij  dr 


(3.16) 


where  T  and  V  are  instantaneous  temperature  and  volume. 
The  virial  sum,  which  is  the  second  term  on  the  right  hand 
side,  can  be  calculated  concurrently  with  the  right  hand 
side  of  (3.1)  since  du/dr  must  be  calculated  anyway. 
Given  the  system  pressure,  we  can  now  calculate  the 
enthalpy  H  from 


H  =  E  +  PV 


(3.17) 


where  V  is  the  system  volume. 

The  radial  distribution  function  g(r)  is  an 
unnormalized  probability  density  function  for  finding  an 
atom  at  a  radius  r  from  another  atom.  It  can  be 
calculated  using 

<N(r  +  Ar/2> 

g(r)  =  - = -  (3.18) 

V(r  +  ^T/2)p 

where  <N(r  +  t\r/2>  is  the  time  average  number  of  atoms 

between  r-^r/2  and  r+^/2,  V(r  +  ^r/2)  is  the  volume  of  the 
spherical  shell  between  these  two  radii,  and p  is  the 
number  density. 

The  mean-square  displacement  MSD(t)  is  the  square  of 
the  average  displacement  of  all  the  atoms  in  the  system  at 
time  t.  For  a  system  containing  N  atoms,  it  is  given  by 


1  _  2 
MSD(t)  =  -  2  [r  (t)  -  r  (0) ] 

N  i  i  i 


(3.19) 


where  r  (t)  is  the  position  of  atom  i  at  time  t. 

1 

The  order  parameter  is  a  measure  of  the  order  of  the 

system.  The  short  range  order  parameter  used  in  this  work 
50 

is  given  by 


S  exp(iK.r  ) 

i»j  ij 


(3.20) 


is  the 


vhere  K  is  a  reciprocal  lattice  vector  and  r 

ij 

vector  between  atoms  i  and  j.  In  this  work  three  values 
of  K  were  used: 


41T 

K  =  —  [1,0,0]  (3.21) 

1  a 

K  =  —  [0,1,0]  (3.22) 

2  a 
Alt 

K  =  —  [0,0,1]  (3.23) 

3  a 

The  order  parameter  results  given  in  Chapter  5  are  the 
averages  of  the  individual  order  parameters  calculated 
with  these  three  reciprocal-lattice  vectors. 

All  the  these  quantities  (except  for  the  total 

energy)  fluctuate  during  the  simulation  and  averages  must 
be  calculated.  For  a  simulation  of  M  time  steps,  the 
average  of  a  property  A(t)  is 


1  M 

<A>  =  -  ^ A(t  ) 

Mm  m 

3.6  The  Parrinello-Rahman  Method 


One  of  the  disadvantages  of  using  conventional 
molecular  dynamics  is  the  that  the  constants  of  motion  are 
the  number  of  atoms  N,  the  volume  V  and  the  total  energy 
E,  which  is  called  a  microcanonical  ensemble.  The  initial 
conditions  (atom  positions  and  velocities)  and  the  system 
volume  specify  N,  V  and  E  at  the  beginning  of  the 


simulation  and  these  values  do  not  fluctuate.  The 

temperature  and  pressure,  however,  fluctuate  with  time 

which  requires  us  to  calculate  a  few  thousand  time  steps 

to  accurately  calculate  their  average  values. 

Consequently,  we  don't  know  the  temperature  and  pressure 

of  the  system  until  after  the  run  is  finished. 

Unfortunately,  experimental  data  for  solids  is  usually 

given  for  a  single  pressure  (i.e.,  atmospheric  pressure), 

which  makes  it  difficult  to  compare  MD  results  with 

experimental  results.  A  way  to  overcome  this  difficulty 

27 

was  proposed  by  Andersen  ,  in  which  one  could  apply  a 

pressure  to  the  system  and  the  volume  would  fluctuate. 

This  was  generalized  by  Parrinello  and  Rahman,  who 

developed  a  method  to  apply  a  general  stress  to  the  system 

28,  51 

by  allowing  the  shape  to  change  .  Parrinello  and 

Rahman's  method  was  implemented  in  this  work  to  allow  us 

to  calculate  thermodynamic  properties  at  a  given  pressure. 

Let  the  system  volume  be  a  parallelopiped,  as  shown 

in  Figure  3-3,  whose  edges  are  the  three  vectors  a,  b,  c. 

Let  r  be  the  position  of  atom  i  in  the  volume.  We  can 
i 

now  define  a  vector  s  such  that 

i 

r  =  h  s  (3.24) 

i  i 

where  h  is  a  3x3  matrix  h  =  {a,  b,  c) .  Since  each  atom  is 

~  k  k  th 

confined  to  the  volume,  0  <  s  <1,  where  s  is  the  k 

~  i  i 

component  of  s  .  We  now  observe  that  we  can  easily  change 


the  volume  and  shape  of  the  system  by  changing  h.  The  key 
to  the  flexible  boundary  method  is  define  equations  of 
motion  for  h  (and  alter  the  equations  of  motions  of  the 
atoms)  such  that  we  can  apply  a  given  hydrostatic  pressure 


to  the  system. 

If  u(r) 

is  a  pairwise 

interatomic 

potential,  these  equations  are 

given  by 

-1 

s  =-m 
i  i 

ZI  (u’/r 

ij 

)(s  -s  )  - 

i  j 

G 

-1. 

Gs  ,  i,j 
“  i 

=  1,2,  ..  .N 

(3.25) 

-1 

h  =  W  (tf 

-  P  )5 

ext 

(3.26) 

where 

o;  =  (h"")  ^ 

(3.27) 

T 

G  =  h  h 

(3.28) 

=  M  h  1 

1 

(3.29) 

T 

V  V  - 

i  i 


H  (uVr 


T 

)r  r 


(3.30) 


atom. 


th 


W  is  physically 


and  m  is  the  mass  of  the  i  atom, 
i 

equivalent  to  a  wall  mass  and  provides  inertia  for  the 

boundaries  of  the  simulation  volume.  Equations  (3.25)  and 

(3.26)  are  valid  when  a  hydrostatic  pressure  is  applied  to 

the  system,  but  they  can  be  generalized  for  an  arbitrary 
51 

external  stress  which  is  not  necessary  in  this  work. 

We  now  make  some  observations  about  the  trajectories 
generated  by  solving  (3.25)  and  (3.26).  The  Hamiltonian 
corresponding  to  equations  (3.25)  and  (3.26)  is 


“U  =  1/2  m  V  +  ^  u(r  )  + 

i  i  i  j>i  ij 

1/2  W  Tr  h^’h  +  pX?  (3.31) 

where  p  is  the  applied  hydrostatic  pressure.  The  first 
two  terms  comprise  the  Hamilitonian  for  conventional 
molecular  dynamics,  i.e,  a  microcanonical  ensemble,  in 
which  the  total  energy  E  is  constant.  The  fourth  term  is 
the  pressure  work  contribution  to  the  enthalpy.  Except 
for  third  term,  the  constant  of  motion  is  then  the 
enthalpy 


H  =  E  +  pfi  (3.32) 

The  third  term  of  (3.31)  is  the  kinetic  energy  of  the 
boundary.  Its  magnitude  can  be  estimated  by  noting  that 
it  contains  nine  degrees  of  freedom  (one  for  each  element 


of  h),  and  the  average  amount  of  energy  in  the  boundary 

will  be  9/2kT,  where  k  is  Boltzmann's  constant.  On  the 

other  hand,  the  first  term  of  (3.31)  is  the  kinetic  energy 

of  the  atoms  and  on  average  is  3N/2kT,  where  N  is  the 

number  of  atoms.  Ignoring  the  third  term  results  in  an 

error  of  3;N,  and  consequently  the  enthalpy  H  is  the 

constant  of  motion  within  this  error. 

One  ambiguity  in  using  the  Parr inello-Rahman  flexible 

boundary  method  is  the  selection  of  the  wall  mass  W.  For 

calculation  of  equilibrium  averages,  the  selection  of  the 

value  of  W  is  immaterial  since  the  averages  are 

independent  of  the  masses,  as  long  as  the  property 

averages  are  accumulated  over  a  few  oscillations  of  the 
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walls.  In  this  work,  we  followed  Andersen  and 

51 

Parrinello  and  and  Rahman  and  used  W  «  20m.  Recently, 

De  Leeuw  et  al  used  a  wall  mass  of  W  =  m  so  that  the 

energy  transfer  between  the  walls  and  atoms  was 

34 

optimized  .  They  also  tested  other  values  of  W  and  found 
that  there  static  results  were  independent  of  the 
selection. 


Molecular  Statics 


4.1  Introduction 

Molecular  statics  (MS)  is  a  method  of  obtaining 
stable  and  metastable  configurations  of  static  systems.  A 
stable  configuration  is  the  lowest  energy  configuration  of 
the  system.  A  configuration  is  metastable  when  it  is  at  a 
local  energy  minimum  such  that  small  perturbations  of  the 
configuration  will  return  to  the  metastable  configuration. 
If  a  metastable  or  stable  configuration  is  used  as  an 
initial  configuration  for  a  molecular  dynamics  (MD)  run, 
the  atoms  would  never  move  provided  they  have  no  initial 
velocities,  i.e.,  the  initial  temperature  of  the  system  is 
zero.  One  can  therefore  consider  MS  as  a  way  of 
calculating  the  properties  of  a  system  at  a  temperature  of 
absolute  zero.  The  two  methods  differ  in  the  desired 
goal,  however.  The  centerpiece  of  molecular  dynamics  is 
the  faithful  solution  of  Newton's  equations  of  motion.  In 
molecular  statics,  one  is  interested  in  determining  stable 
and  metastable  configurations  at  T  =  OK;  the  main  concern 
is  getting  to  those  configurations  as  quickly  as  possible 
and  insuring  the  configuration  is  indeed  stable  or 
metastable. 
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Since  the  earliest  computational  defect  studies  by 
20  18 

Tewordt  ^  al  and  Girifalco  et  al  ,  MS  has  been  the 
preferred  method  compared  to  MD  for  calculating  defect 
formation,  migration,  and  binding  energies.  There  are  two 
reasons  for  this.  The  first  is  that  MS  is  computationally 
much  faster  than  MD.  In  MS  one  is  interested  in  a  single 
configuration  —  the  lowest  energy  configuration  for  a 
given  system.  Once  this  configuration  is  reached,  usually 
after  no  more  than  a  few  hundred  "time"  steps,  the 
computer  run  is  finished.  On  the  other  hand,  using  MD  one 
is  simulating  a  material  at  a  given  temperature  and  one 
must  wait  for  events  to  happen.  For  example,  if  a  vacancy 
jumps  once  a  picosecond  on  the  average,  one  must  calculate 
a  thousand  time  steps  between  each  jump  using  a  time  step 
size  of  0.001  picoseconds.  To  accumulate  good  statistics 
for  calculating  the  average  jump  time  and  other 
properties,  a  hundred  and  preferrably  more  jumps  are 
necessary  which  will  require  at  least  a  hundred  thousand 
time  steps.  In  this  example,  MS  is  three  orders  of 
magnitude  faster  than  MD.  It  also  should  be  pointed  out 
that  the  calculation  of  defect  energies  using  MD  requires 
results  at  different  temperatures,  so  at  least  two  such 
long  runs  must  be  made.  The  second  reason  MS  is  preferred 
over  MD  is  that  data  analysis  is  considerably  easier.  The 
result  of  an  MS  calculation  is  a  single  configuration  of 
atoms  in  which  the  atoms  are  motionless.  This  final 


configuration  can  be  analyzed  on-line,  i.e.,  just  after 
the  run,  and/or  saved  for  future  reference.  In  MD, 
however,  one  has  generated  thousands  of  configurations 
during  a  run.  To  make  matters  worse,  the  atoms  are  moving 
and  are  almost  certainly  displaced  from  their  average 
positions.  This  makes  it  difficult  to  even  locate  the 
defect,  especially  at  high  temperatures. 

Regardless  of  the  advantages  of  MS  over  MD  for  the 

purpose  of  calcuating  defect  energies,  it  must  be  pointed 

out  that  MS  results  are  strictly  valid  only  at  T  =  OK. 

Some  defect  properties  (see,  for  example,  grain-boundary 

50 

melting  transition  work  of  Ciccotti  et  ^  )  cannot  be 

simulated  using  MS.  For  those  properties  which  can  be 

calculated  using  MS,  however,  one  can  argue  that  the 

temperature  dependence  of  the  property  is  overshadowed  by 

the  uncertainty  in  the  interatomic  potential.  If  the 

given  interatomic  potential  will  only  reproduce  the 

migration  energy  of  a  vacancy  within  ten  percent  of  the 

experimental  value,  it  could  be  pointless  to  attempt  MD 

calculations  since  the  variation  of  the  migration  energy 

with  temperature  may  be  less  than  this.  The  only 

extensive  comparison  of  MD  and  MS  defect  calculation  was 
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done  by  Guinan  et  al  .  Guinan  et  ^  calculated  the 
migration  energy  of  interstitials  in  tungsten  using  both 
MD  and  MS,  and  found  the  results  agreed  within  the 
uncertainties  of  the  methods. 


Since  MS  and  MD  are  closely  related,  much  of  the 
discussion  in  the  previous  Chapter  will  apply  to  MS, 
especially  techniques  for  speeding  up  the  calculation  of 
the  potential  function  between  atom  pairs  discussed  in 
Section  3.4.  In  this  Chapter  the  methods  for  getting  the 
metastable  and  stable  configurations  which  will  be 
discussed  since  these  are  peculiar  to  molecular  statics. 
The  following  section  will  cover  previously  used  methods; 
after  that  a  new  method  which  is  computationally  orders  of 
magnitude  faster  than  previous  methods  will  be  presented 
and  all  the  methods  will  be  compared  using  some  simple 
test  cases.  Since  MS  is  used  exclusively  in  this  work  to 
calculate  defect  formation,  binding,  and  migration 
energies,  the  techniques  necessary  to  calculate  these 
quantities  will  be  discussed  in  Section  4. 


4.2  Review  of  Molecular  Statics  Methods 


Different  computational  methods  can  be  used  to 
determine  the  minimum  energy  configuration  of  a  system  of 
atoms.  Before  discussing  these  methods,  however,  it  is 
useful  to  first  point  out  the  characteristics  of  the 
minimization  problem  itself.  The  primary  characteristic 
is  that  the  energy  depends  on  a  large  number  of  variables, 
in  fact  3N  variables  where  N  is  the  number  of  atoms  in  the 
system  and  the  variables  are  the  atom  coordinates. 
Another  characteristic  is  that  we  can  normally  expect  that 
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the  gradient  vector,  which  is  the  first  derivatives  of  the 
the  energy  with  respect  to  the  variables,  is  available 
since  the  derivatives  are  used  to  calculate  the 
interatomic  forces.  On  the  other  hand,  the  second 
derivatives  of  the  energy  with  respect  to  the  variables, 
mathematically  called  the  hessian  matrix  and  proportional 
to  the  force  constant  matrix,  is  not  so  easy  to  calculate. 
The  hessian  matrix  is  not  required  for  MD  calculations  as 
discussed  in  the  previous  Chapter.  Finally,  we  note  the 
system  variables  are  not  bound  by  any  mathematical 
constraints.  Although  the  periodic  boundaries  discussed 
in  Section  3.3  appear  to  constrain  the  atoms  to  the 
simulation  cell,  physically  an  atom  can  leave  the  cell 
because  one  of  its  images  will  enter  from  the  opposite 
boundary.  A  constraint  in  the  mathematical  sense  is 
physically  equivalent  to  the  wall  of  an  infinite  well,  and 
the  periodic  boundaries  allow  us  to  avoid  using  this  type 
of  boundary. 

Historically,  molecular  statics  methods  have  been 
divided  into  two  categories.  The  first  are  methods  in 
which  Newton's  equations  of  motion  are  solved  for  the 
system  of  atoms,  modified  such  that  energy  is  removed  from 
the  system  as  it  approaches  the  minimum  energy  state. 
These  methods  can  be  implemented  by  converting  a  standard 
molecular  dynamics  code  since  the  code  already  solves  the 
equations  of  motion  and  will  be  called  "quasidynamic" 
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methods.  The  other  category  whose  methods  completely 

dispense  with  solving  Newton's  equations  and  use  general 

mathematical  minimization  techniques  to  minimize  the 
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energy.  These  are  often  called  "variational"  methods 

These  methods  converge  more  quickly  to  the  minimum  energy, 

although  some  believe  that  the  first  category  of 

techniques  will  converge  more  reliably  to  the  global 

53 

energy  minimum  and  not  to  some  local  minimum 

Mathematically,  however,  there  is  no  evidence  that  methods 

in  the  first  category  are  more  reliable  those  in  the 

second,  and  this  author  believes  that  in  general  one 

should  used  the  fastest  method  available. 

There  are  two  MS  quasidynamic  methods  which  are  still 

commonly  used.  (These  methods  not  universally  considered 

MS  methods,  but  are  nevertheless  included  here  since  they 

are  used  for  the  same  purposes  as  the  methods  discussed 

below.)  The  principle  behind  the  first  method  is  to  apply 

a  frictional  or  braking  force  to  the  atoms  which  is 
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proportional  to  their  velocities  .  This  will  be  called 
the  "frictional"  method  in  this  work.  The  justification 
for  this  method  comes  from  the  fact  that  the  total  energy, 
the  sum  of  the  kinetic  and  potential  energies,  for  the 
system  is  constant.  When  the  potential  energy  of  the 
system  is  lowest,  the  kinetic  energy  and  hence  the  atom 
velocities  are  highest  and  the  braking  force  slows  the 
system  down  near  this  configuration.  The  second 


quasidynamic  method  involves  setting  the  velocity  of  an 

atom  to  zero  whenever  the  dot  product  of  the  velocity  and 
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acceleration  of  an  atom  is  negative  .  This  will  be 
called  the  "v.a"  method  here.  This  technique  prevents  an 
atom  from  climbing  out  of  a  local  energy  veil.  The  main 
advantage  of  these  methods  is  that  they  can  be  implemented 
by  making  relatively  trivial  modifications  to  an  existing 
MD  computer  code. 

Two  mathematical  methods  which  used  to  minimize  a 
function  of  many  variables  are  the  steepest-descent  and 
conjugate-gradient  methods.  These  methods,  like  the 
quasidynamic  methods  discussed  above,  only  require  the 
gradient  vector  to  be  calculated.  When  one  has  a  working 
MD  computer  code,  the  main  disadvantage  of  the 
steepest-descent  and  conjugate  gradient  methods  is  that 
one  must  make  more  extensive  modifications  to  an  existing 
MD  code  to  implement  them. 

Probably  the  fastest  minimization  schemes  for 
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unconstrained  problems  are  Newton-like  methods  .  (Note: 

Newton-like  mathematical  minimization  methods  are  not 

related  to  Newton's  equations  motion.)  These  techniques 

use  calculated  or  approximated  first  and  second 

derivatives,  i.e.  the  gradient  and  the  hessian,  to  perform 

the  minimization.  Let  r  be  the  atom  positions  at  the  end 

m 

of  the  m  iteration,  and  g  and  [Q  ]  be  the  corresponding 
th  mm 

gradient  and  hessian  matrix.  The  next  iteration  r  is 

m+1 


then  given  by 


-1 


The  Newton  method  has  recently  be  used  by  Wolf  in  grain 
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boundary  statics  calculations  .  Wolf  found  that  only  a 

few  iterations  were  required  to  minimize  grain  boundary 

configuration  which  would  required  hundreds  of  iterations 
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using  gradient  methods 

There  are  some  important  drawbacks  to  this  Newton's 
method,  however.  The  first  is  a  drawback  to  the 
programmer:  it  is  necessary  to  calculate  to  hessian  matrix 
[Q].  This  is  not  a  trivial  task,  but  if  it  only  has  to  be 
programmed  once  and  a  great  improvement  in  computer 
efficiency  can  be  achieved,  then  it  is  probably  warranted. 
The  second  drawback  is  the  requirements  on  computer 


memory. 


[Q]  is  a  3Nx3n  symmetric  matrix.  Assuming  four 


bytes  per  real  number  and  N=500,  4.5  Megabytes  of  computer 
memory  are  required  if  one  takes  advantage  of  symmetry. 
This  amount  of  memory  was  not  available  for  this  work,  but 
most  modern  computers  (including  microcomputers)  can  be  so 
equipped  and  this  alone  does  not  make  the  method 
intractable.  The  third  drawback  is  the  most  important. 
In  molecular  dynamics  and  many  molecular  statics  methods, 
most  of  the  computer  time  is  spent  evaluating  the 
interatomics  forces,  i.e,,  the  function  evaluation.  In 
these  situations  the  number  of  function  evaluations 
required  is  a  good  comparison  point  since  it  is 


approximately  proportional  to  the  total  computer  time. 

This  is  not  true  when  using  Newton's  method  on  large 

systems,  however,  because  the  method  introduces  overhead 

which  is  comparable  to  the  time  required  to  evaluate  the 

forces.  Part  of  the  overhead  stems  from  the  calculation 

of  the  hessian  matrix  [Ql.  The  interaction  of  two  atoms 

generates  a  3-component  force  vector  and  a  3x3  matrix  of 

second  derivatives.  If  the  evaluation  of  the  potential 

has  been  optimized,  the  additional  time  required  to 

construct  the  hessian  is  significant  (recall  that  is  also 

necessary  to  calculate  the  second  derivative  of  the 

potential).  The  other  source  of  the  overhead  is  the 

solution  of  the  set  of  linear  equations  for  r  ([Q]  is 

m+1 

not  actually  inverted).  The  time  required  to  solve  these 

3 

equations  is  proportional  to  N  and  hence  is  important  for 
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large  N 

4.3  A  New  Method 

A  special  implementation  of  a  method  due  to  Fletcher 
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and  Powell  allows  us  to  avoid  the  drawbacks  of  Newton's 
method.  What  is  avoided  in  this  implementation  is:  1) 
storage  of  a  3Nx3N  matrix,  2)  direct  calculation  of  the 
hessian,  and  3)  solution  of  a  set  of  3N  linear  equations. 
We  retain,  however,  the  excellent  convergence  properties 
of  Newton-like  methods.  The  standard  implementation  of 
Fletcher's  method  is  shown  in  the  algorithm  in  Figure  4-1. 


The  key  to  the  method  is  the  calculation  of  a  matrix  H 


which  converges  to  [Q]  ,  the  inverse  of  the  hessian. 

Instead  of  calculating  [Q]  directly  and  solving  a  set  of 

equations  for  r  ,  we  update  H  and  do  a  matrix 

m+1  m 

multiplication. 

The  algorithm  in  Figure  4-1  overcomes  two  of  the 

three  drawbacks  of  the  Newton  method  because  the  hessian 

does  not  have  to  be  directly  calculated  and  a  system  3N 

linear  equations  does  not  have  to  be  solved. 

Unfortunately,  we  still  have  to  store  a  3Nx3N/2  matrix 

which  was  beyond  the  capabilities  of  the  computers  used 

for  this  work.  The  solution  to  this  problem  is  found  by 

noting  that  although  A  and  B  are  3Nx3N  matrices,  all  of 

m  m 

the  information  they  contain  is  also  contained  in  the 
3N-vectors  d  and  z  and  two  scalars  a  and  b  such  that 


m  m 


m  m 


T 

d  d 
m  m 

A  =  -  (4.2) 

m  a 

m 


T 

z  z 
m  m 

B  =  -  (4.3) 

m  b 

m 

where 


^*,*^*»  4*"  f”  '  •'*  •***  •"*  -**  > 
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Fiqure  4-1:  Conventional  Algorithm  for 
~  Fletcher’s  Method 


The  portions  of  the  algorithm  which  are 
replaced  in  the  modified  algorithm  in  Figure  4-2 
are  enclosed  in  the  hashed  rectangles  (from 
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Acton  ) 


I  =  0 

Xo  is  given  a  priori,  hence /o.  go 
Ho  =  I  (the  unit  matrix) 


Sj  =  -Hig,  (a  direction,  always  downhill,  and  approximately 
the  distance  to  the  minimum  along  that  direc* 
_ _ Uqn)^ _ 


C :  Execute  a  subroutine  that  minimizes  /  in  the  direction  s,  -  that 
is,  find  the  (positive)  value  of  k  that  minimizes  f(x^  +  ^,)  and 
call  that  critical  value  of  k  by  the  name  a,.  The  nature  of  this 
subroutine  is  irrelevant  to  the  principal  iteration,  although  for 
efficient  operation  a  somewhat  sophisticated  procedure  is 
desirable.  We  have  given  one  earlier. 

C:d,  saA  (the  incremental  distance  we  will 

move  along  s,) 

C;x,+  ,  =  Xj  +  d,  (our  new  position  near  the  mini¬ 

mum  along  Sj) 

C:/(x,+,),gj+,  (note  that  d|’’  gj+,  =  0,  except  for 

rounding  errors) 

C;yj  =  g(+,-g(  _ _ 


yi 


C.B,  = 


^'•^1*1  ^  H,  +  A,  —  B, 


(the  nu'merafo?  is  the  backward,  or  "J 
matrix,  product  of  two  vectors;  | 
the  denominator  is  an  ordinary  dot  i 
product)  j 

(again,  the  numerator  is  the  matrix  | 
product  of  two  vectors;  the  de-  i 
nominator  is  a  quadratic  form,  a  ' 
scalar)  , 

(it  can  be  proved  that  H^,  having  i 
started  positive  definite,  will  always  ' 
remain  g<»niy^definite)  , 


Is  i  <  n? 


Exit,  too  many  iterations! 


/■+  I«—  Test  for  smallness  of  dj  and  Sj  - Stop  (at  minimum) 

no _  small ! 


b  =  y  H  y 
ni  in  m  m 


^4.6) 


We  can  rewrite  H  as 

m 


m-1 

H  =  I  +  H  (A  -  B  ) 

m  j=0  j  j 


(4.7) 


Two  key  observations  can  be  made  at  this  point.  The  first 

is  that  we  do  not  have  to  store  H  directly,  but  instead  we 

can  store  the  sequence  of  vectors  and  scalars  d  ,  z  ,  a  , 

mm  m 

and  b  ,  since  they  contain  all  of  the  information  in  A 
m  m 

and  B  .  As  long  as  m  <  3N/4,  we  will  be  using  less 

m 

storage  than  the  conventional  implementation.  It  would 

appear  that  the  price  of  the  reduced  memory  requirement  is 

a  sacrifice  in  efficiency  since  to  multiply  the  3Nx3N 

matrix  H  times  a  vector,  we  have  to  actually  multiply  the 
m 

vector  by  2m  3Nx3N  matrices.  This  leads  to  the  second 
observation,  where  we  rewrite  the  multiplication  of  the 
vector  w  by  A  as 


=  {d  {  Id  w}/a  }  } 
m  m  m 


(4.9) 


where  the  braces  indicate  the  order  of  evaluation.  Noting 


that  d  w  is  a  scalar,  we  have  reduced  the  matrix 
m 

multiplication  of  w  by  a  3Nx3N  matrix  to  the  evaluation  of 
the  dot  product  of  two  3N  vectors  and  the  multiplication 
of  a  3N  vector  by  a  scalar.  The  net  result  is  that  this 
implementation  is  actually  faster  than  the  conventional 
method  for  m  <  3N/2,  which  is  definitely  the  case  for  the 
molecular  statics  problems  in  this  work. 

In  Figure  4-2  the  revised  algorithm  is  presented.  It 
was  found  that  this  method  converged  in  about  ten  to  a 
hundred  iterations,  such  that  m  «  3N/4  for  realistic 


4.4  Comparisons  of  Methods 


Some  of  the  MS  methods  discussed  in  the  previous 
sections  were  implemented  to  compare  their  ability  to 
determine  static  configurations.  The  frictional  and  v.a 
methods  were  implemented  since  many  workers  are  familiar 
with  them.  The  steepest-descent  method  was  also 


-77- 


Fiqure  4-2;  Modified  Algorithm  for  Fletcher's  Method 


1=0 

Xo  is  given  a  priori,  hence /g,  go 
Hq  =  I  (the  unit  matrix) 


U  U  ^ 

i  pQ  a  3 


T 

z  z  g 
j-0  b 

j 


C :  Execute  a  subroutine  that  minimizes  /  in  the  direction  s,  -  that 
is,  find  the  (positive)  value  of  X  that  minimizes  /(x,  +  ^i)  and 
call  that  critical  value  of  A  by  the  name  a,.  The  nature  of  this 
subroutine  is  irrelevant  to  the  principal  iteration,  although  for 
eOicient  operation  a  somewhat  sophisticated  procedure  is 
desirable.  We  have  given  one  earlier. 


C  id|  = 

C:x,n  =  X,  +d, 

C  •/(*(♦!)•  Sit  t 
C:y,  “  giti  -  gi 


i  I  j-O  a  b 

i  i 


C:a  •  d  y 
i  i  i 


(the  incremental  distance  we  will 
move  along  s,) 

(our  new  position  near  the  mini¬ 
mum  along  s,) 

(note  that  d,''  ■  g,  4.  |  =  0,  except  for 
rounding  errors) 


T  T 

d  d  y  a  z  y 

j  j  i  W  j  j  I 


Cib  -  y  z 
i  i  I 


Is  I  <  n?  >  Exit,  too  many  iterations! 
no 


/  =  ;  +  1  ^  Test  for  smallness  of  d,  and  s,  - jn*  Stop  (at  minimum) 
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implemented;  the  conjugate  gradient  method  did  not  appear 
to  converge  significantly  faster  in  some  tests,  although 
those  results  are  not  given  here.  Because  of  memory 
limitations,  Newton's  method  could  not  be  implemented  for 
this  test.  Finally,  the  Fletcher-Powell  method  as 
discussed  in  the  previous  section  was  implemented. 

Three  test  cases  were  selected  to  compare  the 
methods.  Case  I  is  the  minimization  of  a  single  vacancy 
system.  A  single  atom  was  removed  from  a  108  atom  fee 
cube,  and  this  was  used  as  the  initial  configuration.  In 
Case  II  the  methods  were  used  to  minimize  a  system 


containing  an  19-vacancy  platelet. 


The  initial 


configuration  was  created  by  removing  19  atoms  forming  a 
hexagon  from  a  (111)  plane  of  a  108  atom  fee  cube.  Since 
more  atoms  are  removed  in  the  Case  II  than  Case  I,  one 
expects  that  more  atoms  will  be  displaced  and  the 
displacements  will  be  larger.  Case  III  is  a  defect 
migration  calculation  for  a  single  vacancy.  In  this  Case 
the  system  is  constrained  above  the  minimum  energy 
configuration  which  makes  minimization  more  difficult  than 
in  Cases  I  and  II  (this  method  is  discussed  in  more  detail 
in  the  following  section).  The  impetus  to  find  a  faster 
MS  technique  was  because  the  migration  calculations  took 
so  long  using  existing  methods. 

The  results  of  the  comparisons  are  shown  in  Figures 
4-3,  4-4  and  4-5.  It  is  evident  that  the  both  the 


steepest  descent  and  Fletcher-Powell  methods  are 
considerably  faster  than  the  frictional  and  v.a  techniques 
in  Figure  4-3  and  4-4.  The  steepest-descent  method  is 


only  about  a  factor  of  two  slower  in  these  two  Cases.  In 
the  migration  calculation  in  Case  III,  it  is  apparent  the 
Fletcher-Powell  method  is  far  and  away  the  best  technique. 


It  converges  an 

order 

of  magnitude 

faster  than 

the 

steepest-descent 

method. 

and  within 

500  steps 

the 

frictional  and  v.a  methods  have  converged  to  about  four 
significant  figures. 

4.5  Defect  Energy  Calculations 

4.5.1  Formation  Energy 

The  defect  quantities  calculated  in  this  work  are  the 
formation,  migration,  and  binding  energies  of  defect 
clusters.  In  this  section  we  discuss  the  use  of  molecular 
statics  to  calculate  these  defect  properties.  The 
calculation  of  the  formation  and  binding  energy  of  a 
defect  requires  the  use  of  molecular  statics  in  a 
straightforward  manner  to  determine  the  metastable 
configurations  and  energies  of  a  defect,  and  then  the 
proper  manipulation  of  the  energies  to  get  the  formation 
and  binding  energies.  To  calculate  the  migration  energy, 
however,  one  must  make  some  minor  modifications  to  a 
standard  molecular  statics  code. 


Figure  4-3:  Comparison  of  MS  Methods  for  Case 
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The  formation  energy  of  a  defect  is  the  difference  in 
energy  between  a  system  S  containing  the  defect  and  an 
otherwise  identical  system  S*  which  does  not.  The 
formation  energy  of  almost  any  defect  can  be  defined  for  a 
system  S*,  such  as  an  interstitial,  vacancy,  defect 
cluster,  etc.  S*  is  normally  selected  to  be  the  ideal  or 
perfect  lattice,  and  the  formation  energy  is  the  energy 
required  to  introduce  a  defect  into  a  perfect  lattice. 
Generally  speaking,  S*  does  not  have  to  be  a  perfect 
lattice;  for  example,  S*  can  be  a  system  containing  a 
grain  boundary  and  the  "defect"  could  be  a  vacancy.  In 
this  instance  the  formation  energy  is  the  energy  required 
to  introduce  a  vacancy  near  a  grain-boundary.  In  the 
context  of  this  work,  S*  will  always  be  the  perfect 
lattice  and  the  defects  will  be  vacancies  and  vacancy 
clusters. 

To  calculate  the  formation  energy  of  a  vacancy  using 
molecular  statics,  it  is  necessary  to  calculate  the 
enthalpies  of  both  a  system  S  containing  a  vacancy  and  the 
corresponding  perfect  lattice  system  S*.  S  and  S*  should 
both  contain  the  same  number  of  atoms  N' .  Let  H  and  H*  be 
the  enthalpies  of  S  and  S*.  The  formation  energy  is  then 

f 

E  =  H  -  H*  (4.10) 

vl 


Unfortunately,  it  is  not  possible  to  calculate  H*  directly 


-8 


using  molecular  statics.  To  simulate  a  perfect  FCC 

lattice,  the  number  of  atoms  N  =  4ijk,  where  i,  j,  and  k 

are  positive  integers,  and  to  simulate  a  lattice  with  a 

vacancy,  N'  =  4ijk  -  1.  In  defect  simulations  one 

normally  uses  i=j=k,  which  further  limits  the  choice  of  N 

and  N' .  The  result  is  that  we  cannot  directly  simulate 

both  S  and  S*  because  there  are  not  corresponding  values 

of  N  and  N' .  The  solution  to  this  difficulty  is  to  take 

advantage  of  the  equivalence  of  each  atom  in  the  perfect 

lattice.  Assuming  the  atoms  are  indistiguishable,  the 

neighborhood  of  each  atom,  i.e.,  the  number  and  positions 

of  surrounding  atoms,  is  identical.  For  this  reason  the 

total  enthalpy  of  a  perfect  lattice  is  proportional  to  the 

number  of  atoms  in  the  system  N.  If  H  is  the  enthalpy  of 

0 

n  perfect  lattice  containing  N  atoms,  and  H*  is  the 
enthalpy  of  a  perfect  lattice  with  N'  atoms,  then 


N' 

H*  =  —  H 
N  0 


(4.11) 


f  N’ 

E  =  H - H 

Vl  NO 


(4.12) 


In  practice,  the  perfect  lattice  system  S  consists  of  N 

0 

atoms  and  the  system  S  containing  the  vacancy  contains 
N'  ®  N-1  atoms.  N  will  be  selected  such  that  N  =  4i-ik. 


S  is  the  system  S  with  a  single  atom  removed,  but  both 

0 

systems  have  the  same  volume. 

It  is  possible  to  rewrite  (4.13)  in  a  form  more 

amenable  for  molecular  statics  calculations.  Let  U,  P, 

and  V  be  the  internal  energy,  pressure,  and  volume  of  S, 

and  U  ,  P  ,  and  V  be  the  corresponding  values  for  S  . 
0  0  0  0 
Then  H  =  U  +  PV  and  H  =  U  +  P  V  .  Assuming  P  =  P  ,  we 

0  0  0  0  0 

get 


f  N-1  r  N-1  1 

E  =U - U  +  P  V - V 

vl  N  0  0  L  N  Oj 

N-1 

=  U  -  -  U  +  <Virial>  (4.14) 

N  0 

where 


N-1  T 

V  - V  (4.15) 

N  oJ 

As  noted  above,  V  =  V  .  The  pressure  is  calculated  from 

0 


<Virial>  =  P 


[ 


P 

0 


Substituting 


z  z: 

6V  i  j/i 
0 

into  (4.15), 


r 

ij 


du(r  ) 

ij 


,  l<i<N 


(4.16) 


<Virial>  = 


du(r  ) 

ij 


i  r  2  r  - 

6N  i  ij  dr 


l<i<N 


(4.17) 


As  noted  above,  every  atom  in  the  perfect  lattice  is 
equivalent.  This  allows  us  to  rewrite  the  sum  over  i  as 


du(r  ) 

ij 


Z  Z  r  - 

i  j/i  ij  dr 


du(r  ) 

Ij 


N  r  - 

j/1  Ij  dr 


,  1  <  i  £  N 


(4.18) 


where  we  have  arbitrarily  chosen  i  =  1  on  the  right  hand 
side.  The  virial  term  then  becomes 


du(r  ) 

Ij 


<virial>  =  -  2-»  r  - 

6  j?<l  Ij  dr 


(4.19) 


The  quantity  <virial>  is  thus  the  virial  contribution  of  a 

single  atom  in  the  perfect  lattice.  For  any  given  lattice 

constant  and  potential  cutoff  radius  we  only  have  tc 

calculate  <virial>  once  since  it  is  independent  of  N. 

Equation  (4.14)  can  be  generalized  to  calculate  the 
f 

formation  energy  E  of  a  cluster  of  n  vacancies  to  get 

Vn 


f  N-n 

E  a  U  -  -  U  +  n  <Virial> 


(4.20) 


where  U  is  the  internal  energy  of  a  system  containing  an 
n 

n-vacancy  cluster. 

Molecular  statics  is  used  to  calculate  the  internal 

energy  values  required  in  Equation  (4.20).  We  normally 

initialize  the  atoms  to  positions  near  the  desired 

configuration  and  use  the  molecular  statics  code  to 

minimize  the  configurational  internal  energy  and  determine 

the  positions  of  the  atoms  in  a  stable  or  metastable 

state.  As  noted  above,  we  only  have  to  calculate  the 

perfect  lattice  energy  U  and  the  virial  term  <virial> 

0 

once  for  a  given  lattice  constant  and  potential,  where 

potentials  are  distinguished  between  different  materials, 

cutoff  radii,  etc.  (Recall  U  scales  linearly  with  the 

0 

number  of  atoms  N,  and  that  <virial>  is  independent  of  N). 
The  internal  energy  U  of  the  lattice  with  the  defect  will 
be  N-dependent  to  some  extent  and  will  obviously  depend  on 
the  choice  of  potential. 

4.5.2  Binding  Energy 


b 


The  binding  energy  E  of  a  cluster  of  n  defects  is 


Vn 

the  difference  in  formation  energies  E 
and  of  n  isolated  defects. 


f 

Vn 


of 


the  cluster 


b  f  f 

E  =  nE  -  E  (4.21) 

Vn  VI  Vn 


Substituting  (4.20)  into  (4.21), 


where  the  virial  term  cancels. 


Molecular  statics  is  used  to  calculate  the  binding 
energy  in  the  same  way  it  is  used  to  calculate  the 
formation  energy,  and  the  same  comments  apply  in  this 
case.  Note  that  it  is  not  necessary  to  calculate  the 
virial  term  <virial>  to  calculate  the  binding  energy. 


4.5.3  Migration  Energy 


In  a  static  lattice,  the  migration  energy  of  a  defect 

is  defined  as  the  energy  barrier  that  the  defect  must 

overcome  to  move  from  one  position  or  configuration  to 

another.  In  this  section  we  will  discuss  a  modification 

to  the  molecular  statics  method  which  allows  one  to 

calculate  the  migration  energies  of  defects. 

Let  us  look  at  a  system  of  N  atoms  with  coordinates 

r  ,  where  i  =  1,,.N,  and  configurational  internal  energy 
i 

U(r  ).  Let  site  A  be  the  current  position  of  a  vacancy, 
i 

and  site  B  be  one  neighboring  sites  to  which  the  vacancy 

A  B 

can  jump.  Coordinates  r  and  r  are  the  atomic 

i  i 


A  A 

coordinates  corresponding  sites  A  and  B,  and  U  =  U(r  ) 

i 


B  B 

and  U  =  U(r  ).  In  Figure  4-6  a  single  vacancy  is  shown 
i 

at  site  A  in  a  two  dimensional  square  lattice.  For  the 


vacancy  to  jump  to  site  B,  atom  1  must  jump  in  the  reverse 
direction  as  shown.  There  are  an  infinite  number  of  ways 
this  jiamp  can  be  accomplished.  No  matter  what  path  is 
taken,  however,  atom  1  must  cross  a  line  drawn 
perpendicular  to  the  jump  direction.  Let  us  define  the 
position  of  this  line  by  a  reaction  coordinate  which  in 
this  simple  case  will  be  the  y-coordinate  of  atom  1  when 
it  is  on  the  line. 


i^(r  )  =  r  (4.23) 

i  i.y 

where  r  is  the  position  of  atom  1  and  r  =  (r  ,  r  ). 

1  1  l,x  i,y 

Every  trajectory  r  (t)  which  goes  from  A  to  B  must  pass 

i 

through  a  given  The  energy  at  the  crossing  is 

U(r  I  r^(r  )  =  i^).  Let  us  define  U(i\^)  as  the  minimum  of 
i  i 

U(r  I  »^(r  )  =  »^) ,  i.e. , 
i  i 

=  min  U(r  I  i^(r  )  =  r^)  (4.24) 

i  i 

U(v\)  is  thus  the  lowest  possible  energy  at  the  crossing 

v^(r  )  =  The  jumping  atom  must  pass  through  all  of  the 

v|^values  between  -1  and  1  to  make  a  successful  jump.  Even 

if  an  atom  follows  the  lowest  energy  path  (defined  when 

U(r  )  =  U(r|(r  ))  ),  it  will  have  to  pass  through  the 

i  '  i 

m 

maximum  of  U(ii)*  We  now  define  the  migration  energy  E 

VI 

by 


(4.25) 


E  =  max  U(/\)  -  U 
VI 

Before  developing  the  method  for  an  fee  lattiee,  we 
need  to  make  some  revisions  in  our  definition  of  the 
reaetion  eoordinate  Effeetively,  serves  as  a 
eonstraint  on  the  system  by  limiting  the  possible 
eonf igurations.  (Note:  this  is  a  eonstraint  in  the 
physieal  sense  beeause  we  are  foreing  the  system  to  be  in 
a  higher  energy  state;  this  is  not  what  is  a  eonsidered  a 
eonstraint  in  the  mathematieal  sense.)  Equations  (4.23) 
and  (4.24)  eonstrain  the  system  by  foreing  atom  1  to  have 
a  speeifie  absolute  y  eoordinate.  In  moleeular  staties 
with  periodie  boundaries,  the  system  ean  be  translated  in 
any  direetion  without  affeeting  the  energy  of  the  system 
beeause  the  energies  are  all  ealeulated  from  relative 
distanees  between  atoms.  If  a  single  atom  is  eonstrained 
to  a  speeifie  absolute  loeation,  the  whole  lattiee  will 
slowly  drift  to  a  minimum  energy  whieh  is  independent  of 
the  eonstraint  (i.e.,  there  is  not  a  one-to-one 
eorrespondenee  between  the  reaetion  eoordinate  and  the 
system  energy).  The  solution  to  this  problem  is  to  define 
a  reaetion  eoordinate  whieh  is  anehored  to  the  atoms  in 
the  lattiee  but  whieh  has  a  similar  physieal  meaning.  A 
reasonable  ehoiee  in  the  two  dimensional  ease  is 


■V  iCTTfTT^T^ 
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r  +  r  +  r  +  r 
2,y  3,y  4,y  5,y 

^(r  )  =  r - (4.26) 

i  l,y  4 

As  long  as  atoms  2  through  5  stay  in  reasonable  positions, 
the  reaction  coordinate  is  well-defined,  is  anchored  to 
the  lattice,  and  has  the  same  physical  meaning  as  Equation 
(4.23) . 

Now  that  we  have  looked  at  the  basic  ideas  of  the 
technique  in  two  dimensions,  let  us  develop  the  method  for 
an  fee  lattice  and  look  at  the  details  of  computer 
implementation.  We  will  first  define  the  reaction 

coordinate  for  the  three  dimensional  case.  In  Figure 
4-7  atom  1  must  move  through  the  "gate"  of  atoms  2  through 
5  to  make  a  successful  jump.  As  before  we  wish  to  define 
a  reaction  coordinate  which  in  this  case  refers  to  a 
plane  instead  of  a  line. We  can  define  a  reaction 
coordinate  to  describe  the  position  of  atom  1  during  the 
jump, 

r  +  r  +  r  +  r 
2  3  4  5 

n(r  )  =  X  u  .  r - (4.27) 

i  1  4 

where  u  =  [1(2/2,  V7/2,  0]  and  ^  is  a  scaling  factor.  We 

-1 

set  =  ('/2a/4)  ,  where  a  is  the  lattice  constant,  so 


that  -1  <  ^  <  1. 

Now  let  us  look  at  the  details  of  the  implementation 
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Figure  4-7:  Migration  Path  for  a  Vacancy  Jump 


of  the  method.  The  basic  problem  is  to  calculate  as 

given  by  Equation  (4.24),  Note  that  an  energy  minimization 

is  required  to  determine  U(ir\^)  from  Equation  (4.24). 

Molecular  statics  is  ideally  suited  for  this,  provided  we 

can  find  a  way  to  enforce  the  constraint  r^(r  )  =  The 

i 

most  natural  way  to  add  the  constraint  is  to  add  an  extra 


term  to  the  Lagrangian  for  the  system, 
unconstrained  system,  the  Lagrangian  L  is 


For  the 


L  =  T  -  U 


(4.28) 


where  T  and  U  are  the  kinetic  and  potential  energies, 
respectively.  Lagrange's  equations  are  then 


d  9L  9L 

dt  dp  dq 

j  j 


(4.29) 


d  dT  -dU 

-  +  -  =0  (4.30) 

dt  3p  aq 

j  j 

where  q  is  a  coordinate  and  p  =  dq  /dt.  The  constrained 

j  j  j 

Lagrangian  LC  is  the  same  as  Equation  (4.28)  with  an  extra 
term. 


LC  =  T  -  (U  +  UC) 


(4.31) 


where  UC  =  UC(i^).  Lagrange's  equations  for  LC  are 


[• 


where  the  first  two  terms  are  the  same  as  in  (4.30),  and 


we  assume  that  UC  is  not  a  function  of 
can  be  expanded 


P  • 

j 


The  third  term 


due 

9q 


due 


£!k 

dq 


(4.33) 


3  3 

Newton's  equations  for  the  constrained  system  are  then 


dv 

km  9  U  due  g  n 

m  =  -  -  —  .  ■■  (4.34) 

dt  3r  d»t 

km  km 

where  q  =r  ,p  =v  ,kis  the  atom  number,  and  m  is 
j  km  j  km 
the  coordinate  (x,  y,  or  z). 

We  now  need  to  specify  ue(f|^).  The  effect  of  ue('^)  on 

the  system  should  be  to  artificially  increase  the 

potential  energy  when  the  system  deviates  from  the  desired 

constrained  configuration.  In  this  case,  UC(»^)  should 

increase  when  f[(r  )  deviates  from  the  desired  value  of  ||. 

i 

Another  feature  of  UC(r|^)  imposed  by  limitations  of 
numerical  integrating  schemes  is  that  it  should  not 
contain  discontinuities.  A  simple  form  for  UC(K|)  which 
has  these  qualities  is 


2 

uc(rL)  =  X  [  ^  ^  J  (4.35) 

where  d  is  the  desired  value  of  f^(r  )  and  {  is  a  constant. 
Figure  4-8  shows  the  effect  of  adding  UC(fl^  defined  by 
(4.35)  to  the  unconstrained  potential  energy  U(rt).  The 


energy  minimum  is  now  near  the  desired  value  of  and  the 
molecular  statics  code  will  converge  to  this 
configuration.  Substituting  into  (4.33)  and  letting  r 


km 


q  . 
j 


sue 


-  =  2  [ri^r  )  -  a] 

9r  i 

km 

and  from  (4.27) 


9r 


(4.36) 


km 


3ri 


km 


(4.37) 


th 

where  j  is  the  m  component  of  j. 
m 

The  simplicity  of  the  method  is  shown  in  Equation 
(4.37).  The  only  atoms  that  contribute  to  the  constraint 
implementation  are  those  specifically  included  in  the 
definition  of  and  only  these  atoms  are  directly 

affected.  Also  note  that  we  can  implement  the  constraint 
technique  independently  of  the  calculation  of  the 
interatomic  forces,  which  means  that  it  can  be 
superimposed  on  those  calculations.  The  molecular  statics 
code  runs  in  basically  the  same  way  as  when  it  is 
unconstrained.  The  difference  is  that  we  add  UC(f^)  from 
Equation  (4.35)  to  the  configurational  internal  (or  total 
potential)  energy,  and  we  add  contributions  to  the  forces 


•  .  "T.  M-w 
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Fiqure  4-8:  Effect  of  Constraint  on  the 
~  Potential  Energy 


U(/|^)  is  the  unconstrained  potential  energy; 
UC(v^)  is  the  additional  energy  added  from  the 
constraint.  Note  hov  the  minimum  of  U(i^)'''UC(^  is 
now  near  the  minimum  of  UC(d) . 


given  by  Equations  (4.36)  and  (4.37).  The  effect  of  these 
additions  is  to  constrain  the  system  near  some  desired 
value  of  »^(r  ) . 


Chapter  5 
Bulk  Properties 


5.1  Introduction 

The  main  thrust  of  this  work  is  the  calculation  of 
the  properties  of  vacancy  defects  in  copper.  As  was 
pointed  out  in  Chapter  1,  many  defect  properties  have  not 
been  experimentally  determined  and  hence  the  computational 
results  are  difficult  to  verify.  The  molecular  statics 
(MS)  and  molecular  dynamics  (MD)  methods  are  powerful 
enough,  however,  to  allow  one  to  calculate  many  different 
physical  properties  using  the  same  interatomic  potential. 
The  purpose  of  the  work  discussed  in  this  Chapter  is  to 
calculate  the  thermodynamic  and  vibrational  properties  of 
the  copper  potential  using  MD  techniques  and  compare  them 
with  experimental  results.  While  correspondence  between 
calculations  and  experiment  in  one  property  does  not 
guarantee  the  same  in  another,  this  way  of  testing  a  given 
potential  is  just  about  the  only  way  of  determining  its 
range  of  applicability. 

To  make  the  test  of  the  Dagens'  potential  even  more 

complete,  similar  calculations  are  also  performed  using 

29 

the  modified  Morse  potential  discussed  in  Section  2.2. 
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By  repeating  the  calculations  with  a  different  potential, 
one  can  better  evaluate  the  sensitivity  of  the  calculated 
property  to  detailed  features  of  the  potential.  As  is 
well-known,  there  are  not  many  materials  for  which  a 
single  potential  is  recommended  for  all  properties,  and  as 
was  discussed  in  Chapter  2,  many  interatomic  potentials 
have  been  used  for  copper. 

5.2  Thermodynamic  Properties 

The  thermodynamic  properties  calculated  in  this 
Chapter  are  the  temperature,  volume,  pressure,  total 
energy  and  enthalpy.  The  behavior  of  these  properties 
with  temperature  allow  us  to  estimate  the  specific  heat 
and  linear  thermal  expansion  coefficient.  To  obtain  some 
information  on  the  structure  of  the  system,  the  order 
parameter  and  mean-square  displacement  are  also 
calculated.  Melting  is  a  first  order  transition,  and 
properties  such  as  the  total  energy  and  enthalpy  are 
discontinuous  at  the  melting  point.  This  transition  can 
be  determined  from  a  series  of  runs,  and  combined  with 
information  given  by  the  order  parameter,  mean-square 
displacement,  and  the  radial  distribution  function,  can  be 
used  to  estimate  the  melting  point  of  the  system.  The 
results  of  these  calculations  are  presented  in  this 
Section. 

The  formulas  and  methodology  for  the  calculation  of 


thermodynamics  properties  were  discussed  in  Chapter  3.  It 
was  noted  that  the  all  of  the  properties  calculated  in 
this  work  are  functions  of  positions  and  velocities.  From 
a  single  configuration  of  the  system,  i.e.,  a  single  set 
of  positions  r  and  velocities  v  of  all  the  atoms,  a  single 
or  instantaneous  value  of  some  property  A(r,  v)  can  be 
calculated.  In  a  constant-volume  or  microcanonical 
simulation,  however,  all  of  the  thermodynamic  properties 
fluctuate  except  for  the  total  energy  E,  the  volume  V,  and 
the  number  of  atoms  N.  If  the  system  remains  in  a  single 
phase,  the  properties  fluctuate  about  an  average  value 
<A>,  where  <A>  is  the  average  of  A(r,  v)  over  all  time. 
As  one  might  expect,  the  value  of  A(r,  v)  between  one  time 
step  and  the  next  are  correlated  (because  of  limitations 
of  the  integrator,  atoms  can  only  move  a  short  distance 
during  a  step).  Consequently,  thousands  of  time  steps  are 
required  to  get  a  reliable  estimate  of  <A>. 

In  Chapter  2  it  was  noted  that  the  structure- 
independent  part  of  Dagens’  noble  metal  potential  is  very 
difficult  to  calculate  and  was  not  available  for  this 
work.  Since  the  simulation  pressure  depends  on  the 
structure- independent  contribution,  only  a  limited  number 
of  derived  thermodynamic  properties  could  be  calculated. 
Furthermore,  it  was  impossible  to  calculate  the  properties 
at  ambient  pressure,  where  experimental  results  are 
usually  quoted.  To  attempt  a  comparison  with  experimental 


data,  thermodynamic  properties  were  calculated  at  three 
specific  volumes  using  conventional  (constant-volume) 


molecular  dynamics.  The  simulation  volumes  V  ,  V  ,  and  V 

Os  1 

correspond,  respectively,  to  the  density  of  copper  at 

ambient  pressure  and  temperature  (1  atmosphere  and  298K); 

the  density  of  solid  copper  at  melting  at  ambient 

pressure;  and  the  density  of  liquid  copper  at  melting 

under  ambient  pressure.  This  is  summarized  in  Table  5-1, 

in  which  T  =  1357K  is  the  melting  temperature  of  copper, 
m 

At  these  volumes  one  can  calculate  the  melting  point  of 

the  potential  and  make  some  tentative  comparisons  with 

real  copper,  assuming  the  thermal  expansion  of  the 

computer  and  real  copper  are  identical.  The 

structure-dependent  part  of  Dagens'  potential,  which  is 

the  part  of  the  potential  used  in  MD  and  MS,  is  also 

volume  dependent.  Even  if  the  structure-independent  part 

of  the  potential  was  available,  strictly  speaking  a 

different  interatomic  potential  must  be  used  at  each 

simulation  volume.  Only  the  single  interatomic  potential 
1 

given  by  Lam  and  specified  by  equation  (2.7)  was 
available  for  this  work  and  was  used  in  all  calculations. 
It  is  not  possible  to  say  how  sensitive  the  form  of  the 
interatomic  potential  is  to  the  total  volume  without 
actually  deriving  the  potential  at  different  volumes. 

The  modified  Morse  potential  used  for  the  parallel 
property  calculations  is  not  volume  dependent  nor  does  it 
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Volume  of 

3 

Unit  Cell  [A  ] 


47.242 


50.489 


53.115 


Corresponding  State  in 

_ Real  Copper 

T  =  298k,  P  =  1  atm 

T  =  T  ,  P  =  1  atm,  solid  phase 
m 

T  =  T  ,  P  =  1  atm,  liquid  phase 
m 


Table  5-1:  Simulation  Volumes  Used  in 
Property  Calculations  with  Dagens' 
Copper  Potential 


have  a  structure- independent  component.  Consequently,  it 
is  possible  to  calculate  all  of  thermodynamic  properties 
of  the  potential  including  the  pressure.  To  take 
advantage  of  this  fact,  the  flexible  border  Parrinello- 
Rahman  method  was  used  for  the  property  calculations.  The 
system  pressure  was  set  to  0  GPa  (effectively,  1 
atmosphere).  Since  experimental  properties  are  normally 
calculated  at  this  pressure,  direct  comparisons  of  all 
properties  could  be  made. 

In  all  of  the  thermodynamic  property  calculations, 
256  atoms  were  used.  In  the  opinion  of  the  author,  this 
is  probably  a  sufficient  number  of  atoms  based  on  previous 
experience  with  calculations  done  for  argon  using  the 
Lennard-Jones  6-12  potential.  The  cutoff  radius  for 
Dagens'  potential  was  set  at  4.9241  A,  which  is  shorter 
that  the  cutoff  radius  used  the  the  defect  calculations 
reported  in  Chapters  6  and  7.  This  cutoff  is  equivalent 
to  1.4a,  where  a  is  the  lattice  constant.  While  the 


author  believes  that  the  results  of  this  work  would  not  be 
significantly  different  it  a  longer  cutoff  radius  was 
used,  this  possibility  was  not  extensively  tested.  As 
noted  above,  thousands  of  time  steps  are  necessary  to 
calculate  reliable  averages  of  the  fluctuating  properties. 
Most  of  the  properties  reported  in  the  following  sections 
were  averaged  over  4000  to  6000  time  steps;  in  some  cases, 
longer  runs  were  made  when  it  was  suspected  that  the 
system  was  in  a  metastable  state.  This  point  will  be 
addressed  in  more  detail  below. 

The  results  of  the  thermodynamic  and  related  property 
calculations  are  presented  in  Tables  5-2  through  5-5.  All 
of  the  thermodynamic  property  values  shown  are  time 
averages  except  for  the  total  energy  E  and  unit  cell 
volume  V  for  Dagens'  potential  and  the  total  enthalpy  H 
for  the  modified  Morse  potential,  which  are  constants  of 
motion  in  these  simulations.  The  numbers  of  time  steps 
over  which  properties  and  their  averages  were  calculated 
are  given  in  Table  5-3  and  5-5  for  each  run.  In  most 
cases  the  averages  were  calculated  over  4000  steps.  Some 
the  properties  are  plotted  in  Figures  5-1  through  5-7  and 
will  be  discussed  below.  The  melting  transition  is 
indicated  in  many  of  the  Figures;  the  discussion  of 
melting  will  follow  the  presentation  of  the  results. 

The  mean-square  displacements  (MSD)  for  the 
pseudopotential  and  the  modified  Morse  potentials  in  the 
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Run  <T>  <P>  <U>  <V>  <E> 

3 

[K]  [KJ/Kg]  [KJ/Kg]  [A  ]  [KJ/Kg] 


1 

0.0 

22.7 

494. 

47.242 

493.94 

2 

280.2 

23.6 

550. 

47.242 

604.95 

3 

285.8 

23.7 

551. 

47.242 

606.91 

4 

921.5 

26.1 

683. 

47.242 

863.77 

5 

1402.2 

28.2 

790. 

47.242 

1065.80 

6 

1904.6 

30.6 

917. 

47.242 

1291.16 

7 

2187.5 

32.6 

1032. 

47.242 

1460.93 

8 

1917.3 

33.5 

1085. 

47.242 

1461.04 

9 

2518.5 

36.5 

1222. 

47.242 

1716.37 

10 

2767.9 

37.8 

1284. 

47.242 

1826.88 

11 

449.9 

15.9 

446. 

50.489 

534.27 

12 

895.6 

17.8 

540. 

50.489 

716.06 

13 

1459.2 

20.5 

672. 

50.489 

958.00 

14 

1576.2 

21.2 

705. 

50.489 

1015.18 

15 

1585.0 

21.3 

704. 

50.489 

1015.20 

16 

1775.1 

22.4 

758. 

50.489 

1106.21 

17 

1827.4 

24.9 

885. 

50.489 

1243.21 

al7 

1758.2 

25.3 

898. 

50.489 

1243.21 

18 

2101.5 

27.2 

982. 

50.489 

1394.43 

19 

2558.1 

29.4 

1079. 

50.489 

1581.44 

20 

773.3 

12.2 

445. 

53.115 

598.77 

21 

931.3 

13.0 

482. 

53.115 

665.11 

22 

1193.4 

14.5 

544. 

53.115 

778.59 

23 

1219.3 

14.7 

551. 

53.115 

790.61 

24 

1294.6 

15.0 

567. 

53.115 

821.33 

25 

1357.1 

18.8 

728. 

53.115 

994.27 

26 

1952.8 

22.1 

869. 

53.115 

1252.35 

27 

2460.3 

24.5 

962. 

53.115 

1444.87 

Table  5-2;  Thermodynamic  Properties  of 
Dagens*  Potential  (Part  I) 


solid  phase  are  plotted  in  Figures  5-1  and  5-2.  The  MSD 
is  calculated  using  equation  (3.19).  It  is  evident  from 
Figure  5-1  and  5-2  that  the  MSD  increases  with  increasing 
temperature.  This  can  be  explained  by  imagining  that  each 


1 

0.0 

1.00 

0.0 

s 

1 

2 

280.2 

0.91 

0.07 

s 

3 

285.8 

0.91 

0.05 

s 

4000 

4 

921.5 

0.75 

0.18 

s 

4000 

5 

1402.2 

0.62 

0.24 

s 

4000 

6 

1904.6 

0.46 

0.44 

s 

4000 

7 

2187.5 

0.27 

1.07 

s->l 

4000 

8 

1917.3 

0.01 

13.97 

1 

6000 

9 

2518.5 

0.00 

9.27 

1 

2000 

10 

2767.9 

0.00 

8.90 

1 

2000 

11 

449.9 

0.88 

0.07 

s 

4000 

12 

895.6 

0.76 

0.14 

s 

4000 

13 

1459.2 

0.57 

0.35 

s 

4000 

14 

1576.2 

0.50 

0.34 

s 

4000 

15 

1585.0 

0.53 

0.38 

s 

4000 

16 

1775.1 

0.45 

0.42 

s 

4000 

17 

1827.4 

0.09 

6.88 

s->l 

4000 

al7 

1758.2 

0.00 

1 

1000 

18 

2101.5 

0.00 

10.67 

1 

3000 

19 

2558.1 

0.00 

16.78 

1 

3000 

20 

773.3 

0.75 

0.16 

s 

4000 

21 

931.3 

0.70 

0.19 

s 

4000 

22 

1193.4 

0.61 

0.32 

s 

4000 

23 

1219.3 

0.62 

0.30 

s 

4000 

24 

1294.6 

0.60 

0.31 

s 

4000 

25 

1357.1 

0.01 

6.78 

1 

4000 

26 

1952.8 

0.00 

10.90 

1 

3000 

27 

2460.3 

0.00 

16.99 

1 

3000 

Table  5-3:  Thermodynamic  Properties  of 
Dagens'  Potential  (Part  II) 


atom  is  inside  of  a  local  energy  well  whose  minimum  is  at 
the  perfect  lattice  site.  As  the  temperature  increases, 


the  atom  is  able  to  climb  further  up  the  sides  of  the  well 
and  hence  the  average  displacement  of  the  atom  also 
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<T> 

<P> 

<U> 

<V> 

<H> 

[K] 

[GPa] 

[KJ/Kg] 

J 

[A  ] 

[KJ/Kg] 

28 

0.0 

0.0 

-1962. 

46.8 

-1962.10 

29 

283.7 

0.0 

-1902. 

47.5 

-1846.18 

30 

295.2 

0.1 

-1900. 

47.5 

-1841.15 

31 

302.1 

0.0 

-1899. 

47.1 

-1838.74 

32 

921.5 

0.1 

-1741. 

49.6 

-1557.00 

33 

1180.6 

0.0 

-1654. 

51.0 

-1417.78 

34 

1329.3 

0.0 

-1600. 

51.9 

-1335.93 

35 

1401.9 

0.0 

-1563. 

52.6 

-1284.82 

36 

1416.0 

0.0 

-1566. 

52.5 

-1284.76 

37 

1226.2 

0.0 

-1444. 

54.9 

-1172.84 

38 

1359.8 

0.0 

-1415. 

55.7 

-1172.78 

39 

1512.5 

0.0 

-1282. 

59.4 

-982.54 

40 

1766.0 

0.1 

-1147. 

65.0 

-797.80 

Table  5-4:  Thermodynamic  Properties  of 
Modified  Morse  Potential  (Part  I) 


increases.  The  MSD  for  Dagens'  potential  shown  in  Figure 

5-1  is  calculated  at  three  different  system  volumes.  Note 

that  as  the  system  volume  increases  from  V  to  V  ,  the  MSD 

0  1 

also  increases  at  constant  temperature.  Dagens'  potential 
is  repulsive  at  these  volumes  (note  that  the  system 
pressure  is  always  positive),  and  hence  the  local 
potential  energy  well  of  an  atom  becomes  shallower  as  the 
volume  increases  which  allows  atoms  to  be  displaced 
further  at  larger  volumes.  It  has  been  estimated  that 
when  the  the  root-mean-square  displacement  exceeds  15-20% 


w  *r  L"*  rt  I 


of  the  nearest  neighbor  distance,  the  solid  phase  is 
unstable  and  a  phase  change  to  the  liquid  phase  should 


- 


''•:a 

'•*  j 


:un  <T>  Order  MSD  Phase  Steps 

2 

[K]  [A  ] 


28  0.0 

1.00 

0.00 

s 

1 

29  283.7 

0.94 

0.06 

s 

4000 

30  295.2 

0.94 

0.06 

s 

3000 

31  302.1 

0.93 

0.07 

s 

2000 

32  921.5 

0.77 

0.27 

s 

3000 

33  1180.6 

0.67 

0.42 

s 

5000 

34  1329.3 

0.60 

0.49 

s 

4000 

35  1401.9 

0.50 

0.75 

s 

1000 

36  1416.0 

0.56 

0.74 

s 

6000 

37  1226.2 

0.01 

10.05 

1 

4000 

38  1359.8 

0.13 

6.58 

s->l 

5000 

39  1512.5 

0.00 

15.85 

1 

3000 

40  1766.0 

0.00 

23.47 

1 

3000 

Table  5-5: 
Modified 

Thermodynamic  Properties 
Morse  Potential  (Part  II) 

of 

60 

occur  .  The 

MSD 

corresponding 

to  the  root 

-mean-square 

displacement  of 

20%  of  the  nearest 

neighbor 

distance 

is 

indicated  in 

Figures  5-1  and  5- 

■2,  and  in  all  cases 

the 

calculated  MSD  exceeded  this  value  without  a  phase  change 
occurring.  It  is  possible  that  the  simulation  run  times 
were  not  sufficiently  long  for  the  phase  change  to  occur, 
although  it  is  noted  that  the  particularly  high  value  of 
0.74  A  in  Figure  5-2  was  averaged  over  7000  time  steps. 
The  problems  associated  with  determining  the  melting  point 
will  be  discussed  in  more  detail  below. 

The  calculated  total  energy  for  Dagens'  potential  is 


izV]  asN 


Figure  5-2;  MSD  vs  Temperature  for  the  Modified 
“  Morse  Potential 


The  dashed  line  indicates  a  root-mean  square 
displacement  of  20%  of  the  nearest  neighbor 
distance. 


plotted  in  Figure  5-3.  The  total  energy  E  is  the  sum  of 


the  kinetic  and  potential  energies  of  the  system  and  is 
calculated  using  equations  (3.11)  through  (3.13).  As 
noted  above,  three  different  system  volumes  were  used  for 
the  thermodynamic  property  calculations  using  Dagens' 
potential.  As  the  volume  of  the  system  is  increased  at 
constant  temperature,  the  total  energy  of  the  system 
decreases.  This  is  caused  by  the  repulsive  nature  of  the 
potential  at  these  densities,  which  causes  the  interaction 
energy  between  atoms  to  decrease  as  they  are  moved  apart 
from  one  another.  Figure  5-3  also  indicates  that  the 
total  energy  increases  with  increasing  temperature  for  a 
fixed  volume.  This  is  caused  not  only  by  the  increase  in 
the  kinetic  energy  of  the  system,  which  is  proportional  to 
the  temperature  as  shown  in  equation  (3.15),  but  also  by 
the  fact  that  the  atoms  on  average  are  further  displaced 
from  the  ideal  lattice  sites  as  the  temperature  increases 
as  shown  in  Figure  5-1  and  5-2.  Consequently,  the  total 
energy  increases  as  the  temperature  increases.  At  some 
point  on  each  curve,  a  vertical  bar  indicates  that  a 
transition  occurs  and  the  total  energy  increases  abruptly. 
This  point  indicates  the  melting  point  at  the  given 
volumes,  which  will  be  discussed  in  more  detail  below.  It 
is  also  noted  that  the  two  data  points,  marked  A  and  B, 
are  slightly  supercooled,  although  the  difference  in 
temperatures  between  these  points  and  their  corresponding 
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melting  points  are  within  the  uncertainties  of  the  melting 
temperatures. 

The  calculated  enthalpy  for  the  modified  Morse 
potential  is  plotted  in  Figure  5-4.  The  enthalpy  is 
calculated  using  equation  (3.17).  All  of  the  modified 
Morse  thermodynamic  calculations  were  performed  at  zero 
pressure.  As  in  Figure  5-3,  the  enthalpy  also  increases 
with  increasing  temperature  for  the  same  reasons. 
Assuming  that  the  indicated  melting  point  of  the  potential 
is  correct,  the  point  marked  C  is  definitely  supercooled. 
This  point  is  Run  37  in  Tables  5-4  and  5-5,  and  the 
properties  were  calculated  over  4000  time  steps. 

The  pressure,  calculated  using  equation  (3.16),  is 
plotted  for  Dagens'  potential  in  Figure  5-5.  It  is 
apparent  from  Figure  5-5  that  a  large  external  pressure  is 
necessary  to  maintain  the  system  at  the  given  volumes. 
The  reason  the  pressure  is  so  large  is  that  the  cohesive 
energy  provided  by  the  electrons  has  not  been  included  in 
the  pressure  calculations  as  discussed  in  Section  3.3. 
Figure  5-5  indicates  that  the  calculated  pressure 
increases  with  increasing  temperature.  This  behavior  is 
caused  by  the  anharmonicity  of  the  interatomic  potential 
as  indicated  in  Figure  2-1.  As  noted  above,  as  the 


temperature  is  increased,  the  displacement  of  the  atoms 
from  their  perfect  lattice  sites  also  increases.  In 
Figure  2-1,  we  note  that  the  potential  energy  increases 


dramatically  with  decreasing  interatomic  radius  compared 


to  the  change  as  the  interatomic  radius  increases.  At 
greater  average  displacements,  the  average  position  of  an 
atom  will  tend  to  be  at  larger  interatomic  radii. 
Consequently,  to  maintain  the  same  average  position  and 
the  same  volume,  an  external  pressure  must  be  applied.  As 
with  the  total  energy,  the  pressure  also  experiences  a 
first  order  transition  at  melting  indicated  by  the 
vertical  bars. 

The  specific  volumes  of  the  modified  Morse  potential 
simulations  versus  temperature  are  plotted  in  Figure  5-6. 
The  behavior  of  the  volume  in  a  constant  pressure 
simulation  is  qualitatively  similar  to  the  behavior  of  the 
pressure  in  a  constant  volume  simulation.  The  volume  of 
the  system  increases  with  increasing  temperature, 
indicating  that  the  potential  is  anharmonic.  The  specific 
volume  also  increases  at  the  melting  point,  indicating 
that  the  system  is  less  dense  as  a  liquid. 

In  Figure  5-7  the  order  parameter  defined  by  equation 
(3.20)  is  plotted  for  Dagens’  potential.  The  order 
parameter  can  be  used  to  determine  the  order  of  a  state 
and  thus  whether  or  not  the  system  is  a  solid  or  liquid. 
The  parameter  is  designed  to  range  between  0  and  1 
inclusive,  such  that  a  1  refers  to  a  perfectly  ordered 
system  (a  cold  fee  lattice,  for  example),  and  0  refers  to 
a  totally  disordered  system.  In  Figure  5-7  the  parameter 

I 


is  plotted  for  the  volume  V  for  the  pseudopotential.  It 

0 

is  interesting  to  note  that  in  the  simulations  at  the 

largest  volume,  V  ,  the  order  parameter  only  decreases  to 

1 

0.6  before  melting,  whereas  it  decreases  to  0.45  at 

volumes  V  and  V  . 

s  0 

Derived  thermodynamic  properties  can  be  determined 
from  the  thermodynamic  data  in  Tables  5-2  and  5-4.  The 
following  properties  are  calculated  as  derivatives  and 
differences  of  the  properties  in  these  Tables: 


-  Specific  heat  at  constant  volume: 


dE 

c  =  — 

V  dT  constant  V 

-  Specific  heat  at  constant  pressure: 


c  =  — 

P  dTl  constant  P 
-  Bulk  modulus: 


dP 

B  =  -V  — 
dV 

-  Bulk  thermal  expansion  coefficient: 


dV 

^  =  — 
T  dT 


(5.1) 


(5.2) 


(5.3) 


(5.4) 


The  latent  heat  of  melting  and  ,  are  given  by  the 

m  m 

differences  in  the  enthalpy  and  the  total  energy. 
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respectively,  of  the  liquid  and  solid  states  at  the 
melting  temperature.  The  volume  change  on  melting  is 
given  by  the  change  in  volume  during  the  melting 
transition,  at  a  constant  temperature.  The  results  of  the 
calculations  of  these  properties  are  shown  in  Table  5-6. 
It  is  noted  that  the  selection  of  the  system  volumes  for 
the  Dagens'  potential  property  calculations  already 
assumes  that  the  thermal  expansion  properties  are  correct, 
and  hence  the  thermal  expansion  coefficient  and  the  volume 
change  on  melting  cannot  be  determined  from  these 
calculations. 

Table  5-6  indicates  that  many  of  the  of  derived 

properties  for  both  potentials  are  in  good  agreement  with 

experimental  results.  The  bulk  modulus  for  Dagens* 

potential  is  within  twenty  percent  of  the  experimental 

value  at  293K.  It  is  noted  that  a  bulk  modulus  of  140  GPa 

for  Dagens*  potential  has  already  been  calculated  at  OK, 

15 

which  gives  even  better  agreement  .  The  specific  heats 
calculated  using  both  potentials  are  within  ten  percent  of 
the  experimental  values.  The  thermal  expansion 

coefficient  for  the  modified  Morse  potential  is  within 
five  percent  of  the  experimental  value  at  293K.  Table  5-6 
therefore  shows  that  both  potentials  are  more  accurate  at 
low  temperatures,  and  near  the  melting  point  the 
calculations  differ  significantly  from  experimental  data. 
Near  melting,  however,  there  are  large  differences  in  the 
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Exp 

142 


Property 

Daqens 

Mod.  Morse 

Exp 

B  [GPa] 

122 

142 

[KJ/Kg] 

m 

74 

246 

205 

c  [J/kg-K] 

V 

393 

367 

c  [J/Kg-K] 

P 

413 

386 

OL  (293K) 

T 

-5 

5.1x10 

4.85x10 

OL  (lOOOK) 

T 

-5 

10.4x10 

6.72x10 

aV 

m 

11% 

5.17% 

T  CK] 
m 

1600 

1480 

1358 

Table  5-6;  Derived  Thermodynamic  Properties 

calculated  and  experimental  results.  The  latent  heat  of 
melting  for  Dagens'  potential  is  less  than  half  of  the 
experimental  value.  The  latent  heat  of  melting  calculated 
using  modified  Morse  potential  is  more  reasonable  and  is 
about  twenty-five  percent  larger  than  the  experimental 
result,  but  there  are  appreciable  discrepancies  in  the 
values  of  the  bulk  thermal  expansion  coefficient  at  lOOOK 
and  the  volume  change  on  melting.  It  is  concluded  that 
the  both  potentials  appear  to  be  fairly  good  at  low 
temperatures,  but  are  not  accurate  at  high  temperatures. 


All  of  the  properties  calculated  in  this  Section  can 
be  used  to  determine  whether  or  not  the  system  is  in  the 
solid  or  liquid  phase  (except  for  the  volume  and  pressure 
in  constant-volume  and  constant-pressure  simulations, 
respectively).  In  plots  of  the  thermodynamic  properties 
such  as  those  in  Figures  5-3  to  5-6,  the  data  will  fall 
into  two  separate  curves  if  both  liquid  and  solid  phases 
were  simulated,  and  if  no  phase  transitions  occured  while 
property  averages  were  calculated.  If  property  averages 
were  calculated  during  the  phase  transition,  the  average 
value  will  fall  between  the  solid  and  liquid  curves.  In 
the  solid  phase,  atoms  vibrate  about  their  perfect  lattice 
positions  and  remain  near  those  positions  as  long  as  the 
temperature  is  not  too  close  to  melting.  Consequently, 
the  mean-square  displacement  (MSD)  will  increase  from  its 
initial  value  of  zero  and  then  fluctuate  about  an  average 
value  which  can  be  uniquely  associated  with  temperature. 
In  a  liquid,  there  are  no  equilibrium  positions  and  atoms 
are  able  to  diffuse  throughout  the  volume,  and  the  MSD 
continuously  increases  (the  increase  is  linear  if  the 
motion  is  diffusive).  Hence,  the  time  evolution  of  the 
MSD  indicates  the  phase  of  the  system.  All  the  properties 
discussed  thus  far  can  only  be  used  to  determine  the  phase 
of  the  system  when  values  of  the  property  are  compared 
between  separate  runs  at  different  temperatures  or  from 
the  behavior  of  the  property  over  time  during  a 


simulation.  A  property  which  can  indicate  the  phase  of 
the  system  from  a  single  instantaneous  value  is  the  order 
parameter.  As  noted  above,  the  order  parameter  defined  by 
equation  (3.20)  has  a  value  of  1  for  a  perfectly  ordered 
system  and  0  if  the  system  is  completely  disordered. 
Although  the  order  parameter  fluctuates  during  a 
simulation,  the  fluctuations  are  small  enough  such  that  a 
single  value  can  often  indicate  the  state  of  the  system. 

An  example  of  how  some  of  these  properties  can  be 
used  to  monitor  a  phase  transition  is  shown  in  Figure  5-8. 
The  internal  energy  U  increases  during  the  run  as  the 
system  leaves  the  low  potential  energy  solid  state  and 
enters  the  higher  energy  liquid  phase.  The  volume  also 
increases  since  the  liquid  is  less  dense  than  the  solid. 
The  order  parameter  decreases  during  the  phase  change; 
note  that  the  fluctuations  in  the  order  parameter  are  much 
less  than  the  fluctuations  in  the  internal  energy  and  the 
volume,  and  that  it  decreases  in  a  smoother  fashion.  The 
MSD  increases  over  the  length  of  the  run  although  it  also 
fluctuates  a  great  deal.  The  advantage  of  having  all  of 
these  properties  available  is  that  they  can  be  used  to 
independently  confirm  one  another.  For  example,  the 
decrease  at  early  times  in  the  order  parameter  is  matched 
by  the  increase  in  the  internal  energy. 

Even  though  one  has  all  of  these  properties  available 
to  determine  the  phase  of  the  system,  the  melting  point  of 


a  potential  is  difficult  to  estimate  accurately  for  a 


number  of  reasons.  One  problem  is  that  the  velocity 

scaling  method  discussed  in  Section  3.3  generates  an 

ensemble  which  has  an  average  temperature  within  about 

five  percent  of  the  desired  temperature.  At  1400K  this 

means  that  the  simulation  temperature  can  be  plus  or  minus 

70K  of  the  desired  temperature.  Hence,  it  is  difficult  to 

accurately  preset  the  temperature  of  the  simulation  near 

the  melting  point  of  the  metals.  Techniques  have  been 

27 

devised  to  overcome  this  problem  ,  but  unfortunately  a 

more  serious  difficulty  exists. 

Ideally,  the  best  way  to  determine  the  melting  point 

of  the  a  system  is  to  calculate  the  free  energy  of  the 

system  in  the  solid  and  liquid  states  at  different 

temperatures.  The  intersection  of  the  two  curves  is  the 

61 

thermodynamic  melting  point  .  Unfortunately,  however, 
the  free  energy  is  very  difficult  to  calculate.  Instead 
of  calculating  the  free  energy,  a  more  commonly  used 
method  is  to  initialize  the  simulation  at  a  given 
temperature  and  phase  and  then  observe  whether  or  not  the 
system  stays  in  that  phase.  If  one  waits  long  enough,  the 
system  will  eventually  enter  the  phase  with  the  lowest 
free  energy.  Unfortunately,  it  may  take  a  long  time  for 
the  transition  to  nucleate,  and  the  information  available 
to  us  from  simulation  cannot  tell  us  how  long  this  will 
take.  (Note  that  the  transition  shown  in  Figure  5-8  took 


over  3000  time  steps).  As  one  moves  away  from  the 
thermodynamic  melting  point,  however,  the  activation 
energy  for  nucleation  of  the  transition  from  a  higher  to  a 
lower  free  energy  phase  decreases  and  the  melting 
transition  can  occur  within  reasonable  simulation  times  of 


a  few  picoseconds.  Consequently,  if  one  continues  to 

raise  the  temperature  of  the  system,  at  some  point  melting 

will  occur.  The  melting  point  under  these  circumstances 

is  called  the  structural  melting  point,  and  the 

temperature  at  which  it  occurs  is  higher  than  the 

61 

thermodynamic  melting  temperature 

Keeping  these  difficulties  in  mind,  the  melting 

points  for  Dagens  potential  for  the  different  simulation 

volumes  are  estimated  to  be  2000K,  1850K,  and  1300K  at  V  , 

0 

V  ,  and  V  ,  respectively.  The  value  of  the  melting  point 
s  1 

corresponding  to  ambient  pressure  is  taken  to  be  the 
average  of  the  two  latter  values,  i.e.,  about  1600K.  The 


5.3  Vibrational  Properties 


The  vibrational  properties  of  a  solid  can  be 

calculated  from  the  positions  and  velocities  of  the  atoms 

at  discrete  times  during  the  simulation.  Although  it  is 

possible  to  calculate  the  vibrational  frequency  spectrum 
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directly  from  the  velocites  ,  the  usual  method  is  to 
calculate  the  velocity  autocorrelation  function  (VACF)  and 
then  take  its  Fourier  transform.  The  VACF  is  given  by 

V  (0)  .  V  (t) 

1  «  i  i 

VACF(t)  =  -  <  Zi  -  >  (5.5) 

N  i  y  (0)  .  y  (0) 

"  i  ""i 

where  v  (t)  is  the  velocity  of  atom  i  at  time  t.  The  <> 
1 

brackets  indicate  a  time  average  which  emphasizes  that  it 

is  necessary  (especially  for  the  number  of  atoms  used  in 

these  simulations)  to  average  over  many  VACF  values 

calculated  for  a  particular  time  t.  As  the  simulation 

0 

proceeds,  the  time  "origins"  denoted  by  t  created 

n 

periodically.  The  atom  velocities  and  their  dot  product 

are  saved  for  each  origin.  To  calculate  the  VACF  at  a 

particular  time  t’,  we  wait  until  that  number  of  time 

0 

steps  has  passed  a  given  origin  n,  i.e.,  until  t-  t  =  t', 

n 

and  calculate  the  numerator  of  (5.5).  In  this  manner  each 
origin  contributes  one  value  to  the  average  for  a  given 
time  t'  of  the  VACF.  For  the  results  that  follow,  100 
origins  separated  by  30  time  steps  were  used,  which  gave 


contributions  to  each  time  t  of  the  VACF.  The  VACF  for 

the  pseudopotential  and  the  modified  Morse  potentials  at 

about  298K  are  shown  in  Figure  5-9.  The  pseudopotential 

VACF  was  calculated  at  V  ;  the  modified  Morse  VACF  at 

0 

ambient  pressure. 

The  Fourier  transform  of  the  VACF  yields  the 

vibrational  frequency  spectrum  of  the  atoms.  The 

vibrational  spectra  for  the  two  potentials  are  plotted  in 

Figure  5-10.  This  plot  also  shows  experimentally  observed 
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spectra  for  copper  at  49K  and  298K  .  The  spectrum  for 

the  pseudopotential  match  the  experimental  observations 

fairly  well,  giving  good  agreement  at  the  positions  of  the 

peaks  and  valleys.  This  confirms  the  conclusions  of 

Upadhyaya  et  who  found  that  the  phonon  spectra  for  the 

pseudopotential  were  in  excellent  agreement  with 
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experimental  data  .  The  spectrum  of  the  modified  Morse 
potential  does  not  match  quite  as  well  and  slightly  misses 
the  high  frequency  peak.  It  is  noted  that  the  system 
volume  using  the  modified  Morse  potential  is  very  close  to 
the  volume  of  real  copper  so  that  this  cannot  be  the 
source  of  the  discrepancy. 


5.4  Discussion 


As  indicated  in  Table  5-6,  both  potentials  appear  to 
perform  better  at  low  temperatures  and  have  difficulty  in 
giving  satisfactory  high  temperature  properties  near 
melting.  The  reason  for  the  excellent  low  temperature 
behavior  is  probably  that  both  potentials,  and  especially 
Dagens'  potential,  reproduce  the  experimentally  calculated 
phonon  spectra  fairly  well,  as  indicated  in  Figure  5-10. 
Of  course,  Dagens*  potential  has  the  additional 
justification  that  it  is  derived  using  the  pseudopotential 
method  and  is  closer  to  a  first  principles  calculation. 
Because  of  this,  it  is  probably  better  to  use  Dagens* 
potential  in  the  low  temperature  regime.  It  is  noted  once 
again,  however,  that  the  Dagens*  interatomic  potential  is 
incomplete  for  these  types  of  calculations  without  the 
structure  independent  contribution. 

Although  the  structure-independent  contribution  to 
Dagens*  potential  was  not  available  for  these 
calculations,  it  is  possible  to  speculate  as  to  whether  or 
not  it  could  make  a  difference.  If  Dagens*  potential  is 
assumed  to  be  correct  at  a  given  volume,  the  volume 
derivative  of  the  structure  independent  energy  is  given  by 
the  difference  in  the  simulation  pressure  and  the 
experimental  pressure  of  copper  under  the  same  conditions. 
This  value,  which  can  be  called  a  pressure  correction  to 


the  simulation  pressure,  should  be  valid  over  all 

temperatures  for  the  given  volume  if  the  potential  is 

correct.  At  V  ,  the  volume  of  real  copper  at  ambient 
0 

temperature  and  pressure,  the  simulation  pressure  is  about 

23.6  GPa  and  hence  the  pressure  correction  at  this  volume 

is  -23.6  GPa.  We  can  calculate  the  pressure  of  real 

copper  at  ambient  temperature  at  the  volume  V  by  using 

s 

the  bulk  modulus  of  real  copper;  this  yields  a  pressure  of 
-9.8  GPa  since  copper  must  be  expanded  to  this  volume. 
The  simulation  at  this  volume  and  temperature  gives  a 
pressure  of  15  GPa;  the  pressure  correction  at  volume  V 

s 

is  thus  -9.8  -  15  =  -25  GPa.  We  can  also  calculate  the 

pressure  correction  at  V  at  the  melting  temperature, 

s 

since  by  definition  of  V  real  copper  is  at  0  GPa.  From 

s 

Figure  5-5  we  find  that  the  simulation  pressure  at  the 
melting  point  of  real  copper,  1357K,  is  about  20  GPa. 
This  gives  a  volume  dependent  contribution  to  the 
potential  of  -20  GPa,  which  differs  from  the  estimate  at 
ambient  temperature  of  -25  GPa  by  20  percent.  Noting  that 
the  bulk  modulus  is  very  sensitive  to  pressures  of  this 
magnitude  since  it  is  calculated  from  a  difference  of  two 
pressures,  it  is  concluded  that  no  volume  dependent 
contribution  could  make  the  pseudopotential  calculations 
match  the  experimental  results. 

The  overall  evaluation  of  the  thermodynamic 
calculations  of  Dagens*  potential  is  that  the  potential 


appears  to  work  well  at  low  temperatures  (less  than 

ambient)  but  is  not  very  successful  at  high  temperatures. 

It  is  interesting  to  note  that  preliminary  thermodynamic 
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calculations  for  Oagens*  silver  potential  ,  which  are 
discussed  in  Chapter  8,  lead  to  the  same  conclusion. 
Although  the  silver  appears  to  be  successful  in  matching 
the  melting  temperature  of  real  silver,  it  does  not  do 
well  with  other  high  temperature  properties.  On  the  other 
hand,  the  low  temperature  properties  calculated  by  the 
silver  potential  are  even  better  than  the  copper  potential 
results. 


3 


Chapter  6 


Small  Vacancy  Cluster  Migration 


6.1  Introduction 


The  primary  factors  which  control  the  concentration 
and  migration  of  defect  clusters  are  their  binding  and 


migration  energies. 


In  general,  the  equilibrium 


concentration  and  the  migration  (or  diffusion)  rate  of  a 
defect  are  described  by  an  Arhennius  relationship  with 
temperature.  For  example,  the  diffusion  coefficient  for  a 
defect  can  be  expressed  as 


D(T)  =  D  exp(-E  /kT) 
0 


(6.1) 


where  k  is  the  Boltzmann  constant,  E  is  the  migration 

energy  for  the  defect,  and  D  is  a  pre-exponential  factor 

0 

which  is  independent  of  or  only  weakly  dependent  on 
temperature.  Since  the  defect  energies  appear  as 
exponents,  they  overshadow  any  minor  temperature 
dependence  of  the  pre-exponential  terms. 

As  discussed  in  Chapter  4,  molecular  statics  can  be 
used  to  calculate  the  formation,  binding  and  migration 
energies  of  defects  in  a  material,  provided  a  suitable 
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interatomic  potential  can  be  found.  The  purpose  of  this 

Chapter  is  to  present  the  results  of  molecular  statics 

calculations  of  defect  energies  using  Dagens'  copper 
15 

potential  discussed  in  Chapter  2.  In  the  following 

section  the  results  of  migration  energy  calculations  of 

the  tri-  and  tetravacanc ies  in  copper  are  presented.  This 

1,  3 

work  extends  the  results  given  by  Lam  et  ^  for  small 

vacancy  defects  in  copper.  A  common  experimental  method 

which  has  been  used  to  determine  small  defect  energies  is 

4 

isochronal  annealing  .  To  make  the  omparisons  with  the 

annealing  experiments  less  ambiguous,  the  experiments  are 

modeled  with  rate  equations  using  the  defect  model  based 

on  the  migration  and  binding  energies  calculated  in  this 

1,  3 

work  and  by  Lam  .  The  annealing  simulations  are 
discussed  in  detail  in  Section  3.  In  the  final  section, 
the  molecular  statics  results  are  compared  to  previous 
computer  calculations  and  experimental  data,  and  the 
applicability  of  the  defect  model  for  copper  is  critically 
discussed. 

6.2  Migration  Energy  Calculations 

As  was  discussed  in  Section  4.5,  the  migration  energy 
between  two  vacancy  defects  which  can  transform  from  one 
to  another  by  a  single  vacancy  jump  is  given  by  the 
magnitude  of  the  energy  barrier  which  must  be  overcome 
during  the  jump.  Many  different  migration  energies  can  be 


-1 


associated  with  a  vacancy  cluster  of  a  given  size  since 

clusters  have  many  different  metastable  configurations. 

Most  of  these  migrations  are  of  no  practical  interest, 

since  a  migration  mechanism  will  be  unlikely  to  occur  if 

it  involves  configurations  which  have  large  formation 

energies  compared  to  other  configurations,  or  if  other 

lover  energy  migration  mechanisms  exist.  Consequently, 

one  is  interested  in  the  minimum  energy  configurations  of 

a  cluster  and  in  the  low  energy  migration  mechanisms  of 

m 

these  configurations.  The  migration  energy  E  of  a 

nV 

n-vacancy  cluster  is  specifically  defined  as  the  energy 
barrier  of  the  lowest  energy  mechanism  which  can  move  the 
the  lowest  energy  configuration  of  the  cluster  from  one 
position  in  a  perfect,  infinite  crystal  to  any  other 
position.  A  perfect,  infinite  crystal  is  specified  to 
eliminate  interactions  with  other  defects  and  surfaces, 
which  are  not  considered  here.  It  is  also  important  to 
specify  that  the  mechanism  be  able  to  move  the  defect 
throughout  the  crystal,  since  some  mechanisms  may  only  be 
able  to  move  the  defect  within  a  limited  region.  The 
strategy  used  to  calculate  the  migration  energy  of  a 
cluster  is  to  calculate  the  binding  energies  of  the 
cluster  configurations  to  determine  the  most  stable 
configuration.  One  then  calculates  the  energy  of  each 
migration  mechanism  which  can  move  the  cluster  from  one 
position  to  another  in  a  crystal. 


All  of  the  defect  energy  calculations  were  done  using 

the  techniques  discussed  in  Section  4.5.  The  system 

consisted  of  a  500  atom  lattice  from  which  the  necessary 

number  of  vacancies  were  removed.  Dagens'  copper 

potential  was  employed  with  a  potential  cutoff  radius  of 

45  -1 

8.5515  A  and  a  Duesbury  factor  of  0.25a  ,  as  discussed 

1,  3 

in  Section  2.3.  Lam  ^  ^  have  already  calculated 

many  of  the  energies  of  small  vacancy  clusters  in  copper 

using  Dagens'  potential,  but  with  larger  systems  (864  to 

2048  atoms)  and  longer  cutoff  radii.  To  insure  that 

system  used  in  this  work  was  sufficiently  large  and  the 

potential  cutoff  radius  was  sufficiently  long,  the  small 

vacancy  defect  energies  already  calculated  by  Lam  et 
1,  3 

al  were  calculated  with  the  system  used  in  this  work. 

Given  the  relatively  large  cutoff  radius  and  system  size, 

and  the  use  of  the  Duesbury  factor,  the  results  of  this 

1,  3 

work  should  correspond  closely  to  those  of  Lam  et  al 

In  Table  6-1  it  is  shown  the  difference  in  the 

1,  3 

calculations  of  this  work  and  those  of  Lam  less  than 
two  percent. 

In  Figure  6-1  a  number  of  copper  tri-  and 
tetravacancy  configurations  are  shown.  Although  many 
configurations  are  possible,  these  configurations  are  of 
particular  interest  because  of  the  role  they  play  in  the 
migration  of  the  clusters  discussed  below.  The  binding 


energies  of  these  configurations  are  shown  in  Table  6-2. 
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Lam  et  al 

This  Work 

m 
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0.81 

IV 

b 

1 
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0.05 

0.05 
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0.55 

0.55 
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0.15 

0.15 
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1 
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0.28 

4V 
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1,  3 

Results 

Lam  et  al 


Cluster 

Binding  Energy  (ev) 

3A 

0.155 

3B 

0.109 

3C 

0.095 

4A 

0.285 

4B 

0.268 

4C 

0.268 

Table  6-2:  Binding  Energies  of  Tri- 
an3  Tetravacancies  in  Copper 

In  Figure  6-2  some  migration  mechanisms  for  the 

trivacancy  are  shown.  In  this  Figure  the  primed  and 

unprimed  states  indicate  the  same  configurations  but 

different  positions  or  orientations.  The  first  migration 

mechanism,  indicated  by  the  3A->3A'  transition,  results  in 

movement  of  the  trivacancy  to  a  different  face  of  an  fee 

tetrahedral  cell.  Although  the  trivacancy  can  reorient 

itself  within  the  tetrahedral  cell,  it  is  not  possible  for 

the  trivacancy  to  escape  the  cell  through  this  mechanism 

and  thus  it  is  not  considered  a  true  migration 
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mechanism  .  The  other  two  mechanisms  in  Figure  6-2, 
3A->3B->3A'  and  3A->3C->3A' ,  do  result  in  migration 
because  through  them  the  trivacancy  is  not  limited  to  a 
small  region  of  the  crystal.  In  these  cases,  the 
trivacancy  moves  through  an  intermediate  metastable 
configuration  before  returning  to  the  stable  configuration 
at  a  new  site.  In  Figure  6-3  configuration  energy  versus 
reaction  coordinate  is  plotted.  The  plot  for  the 


reorientation  mechanism  is  interesting  because  it 


indicates  that  there  is  an  intermediate  configuration 


analagous  to  the  stable  trivacancy  for  the  calculated  for 

3 

silver  by  Lam  ^  al  .  The  plots  for  the  two  true 

migration  mechanisms  yield  no  surprises  and  show  that  the 

migration  which  passes  through  configuration  3B  is 

m 

energetically  favored  and  E  =  0.56  ev,  which  is  about 

3V 

the  same  as  the  divacancy  migration  energy. 

The  migration  mechanism  for  the  copper  tetravacancy 

is  shown  in  Figure  6-4.  The  tetravacancy  passes  through 

three  intermediate  configurations  to  migrate,  which  can  be 

written  4A->4B->4C->4B' ->4A' .  The  reaction  coodinate 

plots  in  Figure  6-5  show  that  the  migration  from 

configuration  4A  to  4B  has  the  highest  energy  barrier,  and 

hence  it  controls  the  energy  for  the  overall  migration, 
m 

Thus,  E  =0.38  ev,  which  is  lower  than  the  calculated 
4v 

migration  energies  of  single,  di-,  and  trivacancies. 


6.3  Annealing  Kinetics 


I 


I 


I 


The  results  of  the  calculations  of  the  defect 

energies  in  copper  are  summarized  in  Table  6-3.  An 

experimental  method  which  has  been  used  to  determine  these 

4  4 

values  in  metals  is  isochronal  annealing  .  As  Balluffi 
points  out,  however,  annealing  experiments  are  not  always 
easy  to  interpret,  especially  when  complicated  processes 
are  involved.  In  copper,  for  example,  only  the  single 
vacancy  migration  energy  is  known  although  the  di vacancy 


i 


X 
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is  believed  to  be  more  mobile  .  Direct  comparison  of 

the  defect  energies  of  large  clusters  with  experimental 

values  is  hence  not  possible  and  another  approach  is 

necessary.  The  tri-  and  tetravacancy  migration  energies 

calculated  in  the  previous  section  complete  the  set  of 

binding  and  migration  energies  of  vacancy  clusters  with 

1,  3 

one  to  four  vacancies  calculated  by  Lam  ^  ^ 

Although  one  cannot  practically  simulate  annealing  using 

molecular  statics  and  dynamics,  the  defect  energies  can  be 

used  in  other  theoretical  and  computational  methods  which 

are  capable  of  doing  modeling  annealing.  In  particular, 
66 

rate  equations  are  very  useful  for  modeling  of  large 

systems  in  which  the  defect  distribution  is  homogeneous 

and  is  often  used  to  interpret  annealing  results  (for 

64 

examples  in  particular  metals,  see  De  Jong  and  Koehler  , 

67  68 

Seeger  and  Frank  ,  and  Collins  et  al  ).  This  gives  one 
the  ability  to  model  defect  annealing  in  copper  and 
compare  the  results  with  experimental  observations.  In 
this  section,  a  general  background  will  first  be  given  on 
annealing  experiments,  followed  by  the  method  and  results 
of  modeling  the  experiments  using  the  defect  energies 
calculated  by  molecular  statics  and  dynamics. 


The  initial  step  of  an  annealing  experiment  is  the 
introduction  of  defects  into  the  material  using,  for 
example,  electron- i rradiat ion  or  quenching.  The  basic 


Energy  (eV) 

Defect  Binding  Migration 


V 

0.82 

1 

* 

V 

0.05 

0.55 

2 

* 

V 

0.15 

0.56 

3 

* 

V 

0.28 

0.38 

4 

Table  6-3:  Binding  and  Migration  Energies 
^  Small  Vacancy  Clusters  in  Copper 


*  1,  3 

(  from  Lam  et  al  ) 


idea  behind  annealing  experiments  is  then  to  activate  as 

few  different  types  of  defects  at  a  time  as  possible  by 

systematically  raising  the  temperature  of  the  specimen. 

Electron-irradiation  creates  single  vacancies  and 

interstitials  in  the  material;  at  low  temperatures  the 

difference  in  the  migration  energies  of  vacancies  and 

interstitials  allows  the  latter  to  migrate  while  the 

former  is  effectively  immobile.  The  general  temperature 

regimes  over  which  each  type  of  defect  anneals  are  called 

69 

Stages.  Within  the  one-interstitial  annealing  model  ,  as 
the  temperature  of  the  material  is  raised,  the 


interstitials  recombine  with  vacancies,  anneal  at  sinks, 
and  form  interstitial  clusters  during  Stage  I.  In  Stage 
II.  interstitial  clusters  coarsen  and  some  interstitial 


reactions  with  impurities  may  occur.  In  Stage  III,  which 

4 

occurs  at  about  280K  in  copper  ,  there  are  no  single 
interstitials  and  the  remaining  single  vacancies  begin  to 
migrate.  Similar  to  interstitials,  the  vacancies  can  form 
vacancy  clusters  and  annihilate  at  the  interstitial 
clusters  formed  during  Stages  I  and  II  and  at  extended 
sinks.  At  higher  temperatures,  two  additional  Stages  are 
found  in  most  metals  in  which  the  vacancy  clusters  coarsen 
and  eventually  all  of  the  original  damage  is  annealed. 

Another  commonly  used  method  of  introducting  defects 
into  a  material  is  quenching.  The  material  is  taken  to  a 
high  temperature,  near  melting,  and  the  relatively  low 
formation  energy  of  vacancies  allows  an  appreciable  number 
of  vacancies  to  exist  in  equilibrium.  The  sample  is  then 
quenched  to  a  very  low  temperature  at  which  the  vacancies 
are  supersaturated  but  also  immobile.  The  temperature  of 
the  specimen  is  then  systematically  increased  as  in 
electron-irradiation  case,  although  of  course  we  do  not 
expect  to  observe  Stages  I  and  II  since  no  interstitials 


are  present, 


The  advantage  of  quenching  over 


electron-irradiation  is  that  only  vacancies  are  formed  in 
quenching  because  the  formation  energy  of  interstitials  is 
too  high  for  them  to  exist  at  an  appreciable  concentration 
even  at  high  temperatures.  On  the  other  hand,  during  the 
quench  it  is  possible  for  vacancies  to  form  clusters, 
which  may  make  it  difficult  to  interpret  the  experimental 
results. 


The  defect  concentration  is  usually  monitored  using 


some  macroscopic  material  property  P  (traditionally  this 

4 

has  often  been  electrical  resistivity  ),  and  analysis  of 

the  changing  concentration  with  time  and  temperature 

enables  one  to  calculate  the  effective  activation  energy 

for  the  annealing  process.  If  annealing  is  controlled  by 

act 

a  single  activation  energy  E  ,  the  annealing  rate  R  is 


given  by 


dP(T)  0  act 

R(t,T)  =  -  =  R  (t)  exp(-E  AT) 

dt 


(6.2) 


where  t  is  time.  R  (t)  depends  on  the  vacancy  defect 

concentrations  at  time  t,  but  is  independent  of 

temperature.  One  method  of  determining  the  activation 

act 

energy  E  from  the  measured  annealing  rate,  called  the 
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change-of-slope  method  ,  is  to  instantaneously  change  the 

sample  temperature.  Let  T  and  T  be  the  sample 

1  2 

temperatures  just  before  and  just  after  an  instantaneous 

temperature  change,  respectively,  and  let  R  and  R  be  the 

1  2 

annealing  rates  just  before  and  after  this  change.  Then 


0  act 

R  exp(-E  /kT  ) 

1  1 

0  act 

R  exp(-E  /kT  ) 


(6.3) 


(6.4) 


Since  the  temperature  change  is  instantaneous,  R  =  R  and 

1  2 

act 

we  can  solve  equations  (6.3)  and  (6.4)  for  E  to  get 

-k  ln(R  /R  ) 
act  1  2 

E  =  -  (6.5) 

1/T  -  1/T 

1  2 

If  a  single  process  (say  single  vacancy  migration)  is 

act 

rate-limiting,  then  a  well-defined  activation  energy  E 
can  be  assigned  to  that  process  (providing  the  process 
itself  can  be  identified).  In  general,  however,  various 
processes  can  take  place  simultaneously;  the  problem  then 
is  identifying  the  rate-limiting  process  or  processes  so 
that  an  effective  activation  energy  can  be  assigned. 

An  annealing  experiment  can  be  modeled  by  solving  a 
set  of  rate  equations  which  describe  the  possible 
reactions  in  the  sample.  Let  Vn  denote  an  n-vacancy 
cluster  where  n  =  1,...,5,  VC  denote  all  larger  vacancy 
clusters,  and  S  denote  sinks.  The  reactions  which  are 
included  in  this  work  are  then  given  in  Table  6-4.  Some 
of  the  possible  cluster-cluster  reactions  have  not  been 
included  in  anticipation  of  the  results  in  which  they  are 
not  significant.  The  sinks  S  in  Table  6-4  can  represent, 
in  the  following  analyses,  both  interstitial  loops  after 
electron-irradiation  and  the  low  concentration  of  sinks 
present  after  quenching.  Although  it  is  possible  for 
appreciable  concentrations  of  large  vacancy  clusters  to 


form,  it  is  not  practically  possible  to  model  all  sizes  of 

clusters.  In  the  model  in  Table  6-4,  all  clusters  with 

six  or  more  vacancies  are  represented  by  VC.  Although  the 

size  of  the  largest  cluster  which  is  explicitly  modeled  in 

Table  6-4  could  have  been  increased  to  a  larger  value,  the 

information  about  these  larger  clusters  is  limited.  The 

rate  constants  for  their  reactions  would  have  to  be 

assumed  without  specific  information  about  the  interaction 

radii  between  the  large  clusters  and  other  defects.  The 

fractions  in  the  VC  reactions  are  an  attempt  to  account 

for  the  reduced  contribution  a  vacancy  makes  to  the 

macroscopic  property  P  when  it  is  a  member  of  a  large 
70 

cluster  .  The  justification  for  these  fractional  values 
will  be  discussed  when  the  macroscopic  property  P  is 
defined  below. 

A  set  of  rate  equations  which  describe  the  reactions 
in  Table  6-4  can  now  be  defined.  Let  the  concentration  of 
a  vacancy  cluster  Vx  be  represented  by  c  ,  and  the  sink 

X 

concentration  by  c  .  The  rate  equations  describing  the 

S 

reactions  in  Table  6-4  are  then 


VI  + 

VI  <=>  V2 

VI 

+  S  ->  S 

V2  + 

VI  <=>  V3 

V2 

+  s  ->  s 

V2  + 

V2  ->  V4 

V3 

+  s  ->  s 

V3  + 

VI  <=>  V4 

V4 

+  s  ->  s 

V4  + 

VI  <=>  V5 

VI  + 

V5  ->  VC 

V2  + 

V5  ->  1.15 

VC 

V3  + 

V5  ->  1.30 

VC 

V4  + 

V5  ->  1.45 

VC 

VI  + 

VC  ->  1.07 

VC 

V2  + 

VC  ->  1.14 

VC 

V3  + 

VC  ->  1.21 

VC 

+ 

> 

VC  ->  1.28 

VC 

Table 

6-4;  Reactions  of 

Vacancy  Defects 

with  one  another 

and 

with  Internal 

Sinks 

2 

=  -2a  c  -  a 

c 

c 

a  c  c  - 

1 

11  1 

12  1 

2 

13  1  3 

a 

c  c  -  a  < 

:  c 

-  a 

c  c  + 

14 

14  51 

1  5 

Cl  1  C 

2b 

c  +  b  c 

- 

2 

2  3  3 

a 

c  c 

SI 

S  1 

c  =  a  c  c  +  1.15  a  c  c  +  1.3a  c  c  + 
C  51  1  4  52  2  5  53  3  5 

1.45a  c  c  +  0.07a  c  c  + 

54  4  5  Cl  1  C 

0.14a  c  c  +  0.14a  c  c  + 

C2  2  C  C3  3  C 

0.28a  c  c 
C4  4  C 


(6.11) 


The  rate  constants  a  are  for  association  reactions 

xy 

between  defects  represented  by  the  concentrations  c  and 

X 

c  .  The  second  subscribt  denotes  the  more  mobile 

y 

component  of  the  two  defects.  The  rate  constants  b  are 

X 

for  the  breakup  or  dissociation  of  the  clusters 
represented  by  concentration  c  into  two  smaller  clusters. 

X 

70 

The  rate  constants  are  selected  following  Johnson  .  They 

consist  of  the  following  three  parts:  an  attempt  frequency 

f  ;  a  geometrical  factor  g  ;  and  an  energy  E  which 
y  X  y 

determines  the  temperature  dependence  of  the  rate 

constant.  The  rate  constants  are  then 


a  =  f  g  exp(-E  /kT) 
xy  y  X  X 

b  =  f  g  exp(-E  /kT) 
X  b,x  b,x  b,x 


(6.12) 


(6.13) 


where  the  subscript  'b'  indicates  a  factor  of  a  breakup 
reaction.  The  values  used  in  this  work  are  shown  in  Table 
6-5.  The  main  uncertainties  are  in  the  geometrical  and 
frequency  factors,  since  these  quantities  are  difficult  to 


experimentally  determine  and  have  not  been  actually 


calculated  here.  A  range  of  values  were  tested  for  these 

70 

factors  and,  confirming  Johnson's  conclusions  ,  we  found 
that  the  results  were  not  very  sensitive  to  them.  The 
energies  are  the  main  quantities  of  interest  here  since 
these  were  calculated  in  the  previous  section  using 
molecular  statics,  and  are  the  values  which  are  normally 
determined  in  annealing  experiments.  The  migration 
energies  are  taken  directly  from  the  molecular  statics 
results.  The  "breakup"  or  dissociation  energy  energy  of  a 
cluster  into  two  smaller  clusters  is  calculated  by  adding 
the  migration  energy  of  the  most  mobile  component  to  the 
difference  in  binding  energies  of  the  original  cluster  and 
the  binding  energies  of  the  components.  For  example,  the 
dissociation  energy  of  a  trivacancy  into  a  single  and 
divacancy,  where  the  divacancy  is  more  mobile  than  the 
single  vacancy,  is  given  by 

brk  b  b  m 

E  =  (E  -  E  )  +E  (6.14) 

3V  3V  2V  2V 

As  pointed  out  above,  the  defect  concentrations  are 
usually  not  directly  measured  in  an  annealing  experiment, 
but  instead  some  macroscopic  property  P  which  is  a 
function  of  the  concentrations  is  measured.  P  has 

traditionally  been  the  electrical  resistivity.  It  is 

believed  that  the  contribution  of  a  small  vacancy  cluster 


Defect  f  miq  brk  miq  brk 

14 

V  1x10  0.82  12  12 

1 

13 

V  2x10  0.55  0.87  84  14 

2 

12 

V  4x10  0.56  0.65  144  14 

3 

12 

V  4x10  0.39  0.69  164  14 

4 

12 

V  4x10  0.76  180  14 

5 

V  180  14 

C 

S  12 


Table  6-5:  Values  of  the  Physical  Parameters 
UseH  in  Equations  (6.12)  and  (6.13) 


to  the  electrical  resistivity  is  approximately 

proportional  to  the  number  of  vacancies  it  contains, 

although  it  is  also  known  the  contribution  per  vacancy 
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decreases  as  the  cluster  grows  .  In  this  model  the 
macroscopic  property  P  is  a  simply  a  linear  combination  of 
the  defect  concentrations  given  by 


-157- 

The  decrease  in  contribution  for  larger  vacancies  is 

accounted  for  by  the  fractional  terms  in  the  V  reactions 

C 

in  Table  6-4.  Although  this  model  is  no  doubt  much 
simpler  than  what  occurs  in  real  metals,  the  model  does 
predict  the  correct  remaining  resistitity  after  annealing 
of  a  quenched  sample. 

It  is  difficult  to  simplify  the  equations  (6.6)  to 

(6.11)  without  introducing  unrealistic  assumptions. 

Because  of  this,  it  was  decided  to  solve  the  equations 

(6.6)  to  (6.11)  nummerically.  The  parameters  which  were 

varied  between  different  runs  were  the  initial  single 

vacancy  concentration  c  (0)  and  the  sink  concentration  c  . 

1  S 

The  ranges  of  these  values  covers  the  spectrum  of  standard 

annealing  experiments. 

The  primary  quantities  we  wish  to  calculate  in  these 

runs  are  the  temperature  T  at  which  the  annealing  rate 

ann 

is  greatest  and  the  average  activation  energy  calculated 

using  the  change-of-slope  method.  Figure  6-6  shows  the 

-5 

results  of  an  annealing  run  with  c  (0)  =  10  and 

1 

-6 

c  =  10  .  In  these  plots  we  show  1)  the  "resistivity"  P 

S 

versus  temperature;  2)  the  temperature  derivative  of  P, 

T*dP/dT;  and  3)  the  activation  energy  calculated  at  the 

temperature  changes.  Note  that  as  the  temperature  is 

increased,  a  point  is  reached  at  which  the  annealing  rate 

is  maximum.  The  temperature  at  this  point  is  T  for 

ann 


••  ••  V  %'  V  •«' 


•  «■  *  ^  , 


this  particular  run.  Since  the  effective  activation 
energy  may  vary  with  temperature  during  the  run,  the 
average  activation  energy  is  calculate  from  the  energy 


values  between  about  95%  and  5%  of  the  initial  resistivity 

R  . 

0 

The  overall  results  of  the  rate  equations  solutions 

are  shown  in  Table  6-6  and  6-7.  Table  6-6  indicates  that 

the  activation  energies  do  not  follow  an  easily 

predictable  pattern.  In  general,  the  effective  activation 

energy  is  fairly  close  to  the  single  vacancy  migration 

energy  of  0.82  eV.  This  indicates  that  single  vacancy 

migration  controls  the  annealing  rate  in  most  runs.  The 

-5 

primary  exception  is  the  case  in  which  c  (0)  =  10  and 

1 

-6 

^5=  10  .  This  run  is  shown  in  Figure  6-6.  In  the 

beginning  of  the  run  the  activation  energy  is  relatively 
low  and  then  increases  to  an  energy  similar  to  the  average 
values  in  other  runs.  This  indicates  that  initially 
another  process  in  addition  to  single  vacancy  migration  is 
also  controlling  the  annealing.  At  early  times  the 
primary  pathways  for  the  annealing  of  vacancies  is  the 
combination  of  single  vacancies  to  form  divacancies,  which 
then  migrate  to  sinks.  At  later  times,  the  single  vacancy 
concentration  decreases  such  that  a  single  vacancy  will 
annihilate  at  a  sink  before  it  can  combine  with  another 
vacancy.  Hence,  at  early  times  in  this  particular  run. 


divacancy  migration  partially  controls  the  annealing  and 
the  divacancy  migration  energy  affects  the  calculated 
effective  activation  energy. 


-log(c  (0)) 

1 

4 

-logic  ) 

S 

5 

6 

8 

4 

OT8(5 

077^ 

— Tm — 

5 

0.81 

0.68 

0.74 

6 

0.81 

0.79 

0.81 

Table  6-6;  Average  Effective  Activation  Energies 

of  Anneals  (eV) 


The  behavior  of  the  maximum  annealing  temperature 

T  with  sink  and  initial  vacancy  concentration,  shown  in 
ann 

Table  6-7,  is  easier  to  explain  than  the  behavior  of  the 

average  annealing  temperature.  We  find  that  at  high  sink 

concentrations,  T  is  relatively  independent  of  initial 

ann 

vacancy  concentration,  but  at  low  sink  concentrations  the 

initial  vacancy  concentration  does  affect  T  .  At  high 

ann 

sink  concentrations,  the  primary  annihilation  mechanism 

for  vacancies  is  absorption  by  sinks.  This  is  a  first 

order  reaction  since  its  rate  equation  term  is 

proportional  to  the  single  vacancy  concentration,  given  by 

a  c  c  .  The  probability  that  an  individual  vacancy  will 
SI  1  S 

interact  per  unit  time  is  simply  the  rate  equation  term 

divided  by  c  ,  i.e.,  a  c  .  Hence,  the  probability  that  a 
1  S  S 

particular  individual  vacancy  will  interact  with  a  sink  is 


independent  of  vacancy  concentration.  As  we  increase  the 


vacancy  concentration,  we  expect  more  vacancies  to 

interact  with  sinks  per  unit  time,  but  the  peak  annealing 

rate  should  occur  at  the  same  temperature.  At  low  sink 

concentrations,  instead  migrating  directly  to  sinks, 

single  vacancies  form  divacancies  first  and  which  either 

migrate  to  a  sink  or  perhaps  form  larger  clusters.  The 

rate-limiting  step  is  the  vacancy-vacancy  interaction 

2 

which  is  a  second  order  reaction  given  by  a  c  . 

11  1 

(Strictly  speaking,  the  overall  reaction  will  not  be  quite 

second  order,  since  the  dissociation  of  divacancies,  a 

first  order  reaction,  will  also  play  a  role).  The 

probability  that  an  indivividual  vacancy  will  interact 

with  another  vacancy  is  proportional  to  a  c  .  Hence,  the 

11  1 

higher  the  concentration,  the  faster  the  reaction  will 
occur  and  T  will  occur  at  a  lower  temperature. 


ann 


-log{c  (0)) 

1 

4 

-log(c  ) 

S 

5 

6 

8 

4 

282 

285 

287 

287 

5 

291 

306 

325 

325 

6 

291 

312 

335 

383 

Table  6-7;  Temperatures  for  Maximum  Annealing  [K] 
6.4  Discussion 


Compared  to  calculations  of  the  binding  energies  of 
vacancy  clusters,  there  have  been  very  few  calculations  of 
their  migration  energies  for  other  than  divacancies.  The 


primary  reason  for  this  is  that,  although  single  and 

divacancy  migrations  in  fee  metals  occur  through  simple 

and  predictable  mechanisms,  the  large  number  of  possible 

configurations  of  large  clusters  makes  enumeration  of 

their  migration  paths  difficult.  As  was  shown  in  Section 

2,  it  is  necessary  to  first  determine  the  stable  cluster 

conf iguration(s) ,  and  then  calculate  migration  energies  of 

favorable  migration  paths.  The  first  work  in  this  area 

was  done  by  Johnson  on  vacancy  clusters  in  nickel  using 

22 

the  Johnson  I  potential  .  After  calculating  the  binding 
energies  various  vacancy  clusters,  Johnson  discovered  that 
the  energies  could  be  accurately  estimated  by  bond 
counting  (creation  of  a  vacancy  breaks  twelve  nearest 
neighbor  bonds).  Bond-counting  suggests  that  the  most 
compact  clusters  (voids)  are  the  most  stable,  which  is 
definitely  not  the  case  for  the  copper  potential  in  this 
work.  As  is  shown  in  the  following  Chapter,  beyond 
clusters  of  ten  or  more  vacancies,  stacking  fault 
tetrahedra  and  hexagonal  loops  are  formed  and  voids  are 
not  observed.  Johnson  also  noted  that  the  calculated 
migration  energies  of  vacancy  clusters  could  be  estimated 
by  the  divacancy  migration  energy  added  to  one-half  the 
difference  in  the  binding  energies  of  the  initial  and 
final  configurations.  Since,  to  migrate,  a  large  cluster 
must  partially  disassociate,  i.e,  move  through  a 
higher-energy  configuration,  the  large  cluster  migration 
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energies  are  always  higher  than  the  divacancy  migration 

energy.  This  is  obviously  not  the  case  in  this  work  in 

which  the  tri-  and  tetravacancy  were  at  least  as  mobile  as 

the  divacancy.  The  Johnson  I  potential  probably  yields 

such  simple  formulas  for  vacancy  properties  because  it 

contains  only  first  nearest  neighbor  contributions. 

Dagens'  potential  used  in  this  work  is  truncated  at  the 

twelfth  nearest  neighbor  distance,  and  the  lattice  must  be 

held  together  by  applying  an  external  stress.  Although 

simple  formulations  have  been  used  to  estimate  the  vacancy 

72 

cluster  migration  energies  ,  the  conclusion  here  is  that 

vacancy  migration  energies  are  very  sensitive  to  the 

potential  and  must  be  calculated  using  molecular  statics. 

The  most  interesting  aspect  of  the  tri-  and 

tetravacancy  migration  results  is  that  the  trivacancy  is 

about  as  mobile  as  divacancy,  and  the  tetravacancy  is  more 

so.  Trivacancies  are  believed  to  be  more  mobile  than 

68 

single  and  divacancies  in  some  metals  such  as  gold  and 
4 

aluminum  based  on  experimental  results,  but  experimental 

studies  have  not  been  able  to  determine  migration  energies 

for  divacancies  or  larger  vacancy  clusters  in  copper. 

There  is  evidence  for  a  vacancy  cluster  in  copper  which  is 

more  mobile  than  the  single  vacancy,  but  the  number  of 

vacancies  in  this  cluster  has  not  been  unambiguously 
5,  6,  7  5 

identified  .  Wienhold  et  al  did  annealing  studies 

of  electron-irradiated  copper  using  doped  and  undoped 
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samples.  Doping  was  used  in  electron-irradiation 

annealing  experiments  to  suppress  the  formation  of  vacancy 

5 

clusters.  Wienhold  ^  al  found  that  the  doped  and 

undoped  samples  both  yielded  the  same  activation  energy  of 

about  0.71  ev,  which  indicated  that  single  vacancies  have 

an  activation  energy  of  about  0.71  ev.  On  the  other  hand, 

the  annealing  of  the  doped  sample  was  governed  by  first 

order  kinetics  and  the  undoped  sample  by  second  order 

kinetics.  First  order  kinetics  in  the  former  sample 

indicates  that  the  single  vacancies  are  annihilating  at 

sinks.  Second  order  kinetics  indicates  that  the 

rate-controlling  step  in  the  annealing  is  the  association 

of  a  defect  with  another  defect  of  the  same  type  (e.g., 

two  single  vacancies  combining  to  form  divacancies). 

Hence,  clusters  must  be  forming,  but  apparently  they  do 

not  affect  the  effective  activation  energy  of  the 

5 

annealing.  As  Wienhold  et  al  point  out,  this  can  happen 

if  the  larger  clusters  (presumably  divacancies)  have  a 

migration  energy  of  about  0.71  ev,  or  if  they  are  so 

mobile  that  they  do  not  limit  the  rate  of  annealing  (i.e., 

they  combine  with  sinks  or  other  vacancies  almost 

instantly  on  formation  because  they  migrate  so  quickly). 

5 

The  latter  possibility  was  selected  by  Wienhold  et  al  as 

the  more  likely  case. 

4 

Balluffi  has  reviewed  the  annealing  studies  for 


A 


io 


various  metals  including  copper.  The  main  conclusion  of 


these  studies  is  that  the  calculated  activation  energy 
using  quenched  materials  is  relatively  constant  and  there 
is  no  indication  of  the  effect  of  divacancies  on  the 
activation  energy. 

The  results  of  the  annealing  runs  summarized  in 

Tables  6-6  and  6-7,  however,  do  indicate  that  the  larger 

clusters  do  affect  the  effective  annealing  activation 

energy.  In  one  particular  case,  it  was  found  that  an 

activation  energy  of  almost  0.2  eV  less  than  the  single 

vacancy  migration  energy  was  calculated.  Although 

experimental  uncertainties  can  mask  some  variation  in  the 

activation  energy,  it  is  difficult  to  believe  that  a  range 

of  0.2  eV  (i.e.,  close  25  percent)  would  be  missed.  It 

should  be  noted  that  the  very  low  sink  concentration  runs 

are  probably  closest  to  quenching  studies  the  activation 

energy  was  fairly  constant.  The  high  sink  concentration 

runs  correspond  to  electron-irradiation  of  doped  samples, 

and  in  these  cases  the  activation  energy  is  also  constant. 

-5 

The  mid-range  sink  concentration  runs,  with  c  =10  and 

S 

-6 

10  ,  probably  correspond  to  electron-irradiation  of 

undoped  copper,  and  it  is  in  these  cases  where  calculated 

and  experimental  results  differ. 

The  results  in  Table  6-7  follow  the  behavior  in  doped 

5 

and  undoped  copper  reported  by  Wienhold  .  At  high  sink 
concentrations,  it  is  found  that  the  primary  reaction  is 


first  order,  but  at  low  concentrations  a  slower  reaction 


is  also  involved.  This  is  not  surprising  since 
divacancies  are  included  in  the  reaction  scheme  in  Table 
6-4,  but  it  does  show  that  neither  larger  vacancies  nor 
the  low  binding  energy  of  the  divacancies  suppress  this 
behavior. 

Other  experimental  evidence  for  highly  mobile  vacancy 

clusters  in  copper  come  from  perturbed-angular-correlation 
6,  7 

(PAC)  studies  .  PAC  is  a  useful  method  because  only 

73 

defects  which  are  mobile  are  identified  .  PAC  studies  on 

copper  indicate  that  two  vacancy-like  defects  are  mobile 
6,  7 

in  stage  III  .  Although  other  mobile  defect  clusters 

were  not  found,  it  is  possible  that  the  highly  mobile  tri- 

and  tetravacancies  calculated  in  this  work  would  be 

invisible  to  PAC.  The  highest  tri-  and  tetravacancy 

concentrations  observed  in  the  annealing  runs  in  this  work 
I  -9  -11 

were  10  and  10  ,  respectively,  which  may  be  too  low  to 

73 

be  detected  .  As  Wichert  points  out,  defects  are  also 

invisible  to  PAC  if  the  defects  do  not  bind  to  the  probe 

atoms,  or  if  the  defect-probe  configuration  has  cubic 
73 

symmetry  .  On  the  other  hand,  based  on  this  work,  we 
cannot  say  whether  or  not  the  mobile  vacancy  cluster 
identified  by  PAC  is  in  fact  the  divacancy. 

The  general  conclusion  of  the  small  defect  cluster 
calculations  is  that  tri-  and  tetravacancies  are  more 

mobile  than  single  vacancies  and  at  least  as  mobile  as 


divacancies  in  copper.  Although  these  results  are 

interesting  in  themselves,  it  is  difficult  to  make 

comparisons  with  experiment  because  the  data  on  tri-  and 

tetravacancies  in  copper  is  nonexistent.  Kinetic  modeling 

of  defect  annealing  in  copper  indicates  that  small  vacancy 

clusters  should  affect  the  calculated  effective  activation 

energy  at  some  defect  and  sink  concentrations,  but  in  most 

cases  the  effective  activation  energy  is  constant  which  is 

4 

supported  by  experiment  .  Since  it  is  possible  to 

question  the  ability  of  the  rate  equations  to  model 

5,  70,  74 

annealing  for  a  number  of  reasons  ,  a  final 


conclusion  cannot  be  made  in  this  respect. 


Large  Clusters 


7.1  Introduction 

The  small  vacancy  defects  investigated  in  Chapter  6 

are  mobile  at  temperatures  as  low  as  280K  in  copper.  If 

vacancies  are  created  at  temperatures  higher  than  about 

280K,  or  if  the  temperature  of  the  specimen  is 

subsequently  raised  above  this  temperature,  the  vacancies 

will  migrate  until  they  combine  with  other  defects  or  are 

75 

absorbed  by  sinks.  Experiments  in  quenched  and 

5 

electron- irradiated  copper  in  Stage  III,  during  which 

vacancies  are  mobile,  indicate  that  much  of  the  initial 

damage  at  the  beginning  of  the  Stage  is  not  annealed, 

which  suggests  that  large  vacancy  clusters  are  formed. 

These  clusters  grow,  not  by  migrating,  but  by  capturing 

single  vacancies  and  small  mobile  vacancy  clusters.  The 

clusters,  if  sufficiently  large,  can  be  observed  by 

8 

transmission  electron-microscopy  ,  and  it  is  found  that 
large  vacancy  clusters  are  in  the  form  of  hexagonal 
dislocation  loops  and  stacking  fault  tetrahedra  (SFT). 

The  purpose  of  this  Chapter  is  to  present  and  discuss 
the  results  of  molecular  statics  calculations  of  the 


-169- 

properties  of  these  large  vacancy  clusters  in  copper. 

76,  77 

Previous  copper  calculations  using  empirical 

potentials  have  determined  the  critical  size  at  which  the 
vacancy  platelets  collapse  into  loops  or  SFT.  However,  no 
calculations  have  been  performed  to  determine  which  of  the 
shapes  are  the  the  most  energetically  stable.  Since  the 
stacking  fault  energy  is  an  important  factor  in 
determining  the  stable  large  cluster  configuration,  it  is 
calculated  first  using  Dagen's  copper  potential.  The 
displacement  fields  and  binding  energies  of  various  shapes 
of  large  clusters  are  then  evaluated.  In  the  final 
section,  these  results  are  compared  with  previous 
calculations  and  experimental  observations. 

7.2  Stacking  Fault  Energy 

It  has  traditionally  been  difficult  to  calculate 

accurate  stacking  energies  using  interatomic  potentials 

and  computer  simulation.  Probably  the  most  important 

reason  for  this  is  that  the  stacking  fault  energy  depends 

on  the  long  range  part  of  the  potential,  i.e,  third 

nearest-neighbor  distances  and  greater,  which  is  not 

accounted  for  very  well  by  empirical  potentials. 

(Empirical  potentials  can,  however,  be  fitted  to  the 

76,  32 

stacking  fault  energy  ).  For  example,  Johnson's  iron 

and  nickel  potentials  only  include  second  nearest-neighbor 

interations,  and  so  have  a  stacking  fault  energy  of 
78 


zero 


The  stacking  fault  calculations  were  performed  by 
setting  up  an  fee  lattice  which  was  rotated  so  that  the 
[111]  direction  was  pointing  in  the  z-direction,  as  shown 
in  Figure  7-1.  This  was  done  because  a  stacking  fault  is 
a  planar  defect,  and  one  can  simulate  an  "infinite"  fault 
by  taking  advantage  of  periodic  boundaries.  The  lattice 
then  consists  of  (111)  or  close-packed  planes,  which  can 
be  stacked  in  any  order.  To  apply  periodic  boundary 
conditions  to  this  rotated  lattice,  the  number  of  planes 
required  to  simulate  a  perfect  lattice  must  be  3n,  where  n 
is  an  integer,  or  3n-l  to  simulate  a  lattice  with  a  single 
stacking  fault.  For  most  of  the  calculations  14  (111) 
planes  were  used,  and  the  boundaries  in  the  x  and  y 
directions  were  made  larger  than  twice  the  cutoff  distance 
for  the  potential. 

The  results  of  the  stacking  fault  calculations  are 

shown  in  Table  7-1.  Also  included  in  this  Table  are 

results  from  calculations  by  Cotterill  and  Doyama  using 

12 

Born-Mayer  and  Morse  potentials  for  copper  .  The 

Born-Mayer  potential  yields  almost  negligible  stacking 
fault  energies.  The  Morse  potential  yields  good  results 
when  shorter  cutoff  radii  are  used,  which  is  probably 
fortuitous  since  the  Morse  potential  is  not  designed  to 
account  for  long  range  interactions.  The  results  for 
Dagens'  potential  are  similar  and  match  the  experimental 
results  at  short  cutoff  radii  of  less  than  12  A.  At 


Figure  2"1*  Schematic  of  System  for 
Calculating  the  Stacking 
Fault  Energy 
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longer  cutoff  radii,  the  stacking  fault  energy  becomes 


negative  which  means  that  the  HCP  lattice  is  more  stable 

15 

than  the  FCC  lattice.  Since  Dagens  has  shown  that  the 

FCC  lattice  is  the  more  stable  than  the  HCP  lattice  for 

the  copper  potential,  the  reason  for  this  discrepancy  is 

1 

most  likely  that  the  fit  presented  by  Lam  et  al  for  the 
interatomic  potential  is  inaccurate  at  the  long  range. 


Potential 

Cutoff 

2 

Energy  [mJ/m  ] 

This  work 

8.5155  A 

65.3 

12.0878  A 

23.3 

14.4194  A 

-22.0 

>14.4194  A 

negative 

12 

Born-Mayer 

176  atoms 

0.007 

4000  atoms 

0.007 

12 

Morse 

176  atoms 

30.8 

200  atoms 

17.9 

4000  atoms 

0.422 

76 

Exp. 

~70 

Table  7-1:  Comparison  of  Stacking  Fault 
Energy  Calculations 


7.3  Large  Cluster  Configurations  and  Energies 


For  consistency,  calculations  of  large  vacancy 

cluster  configurations  and  energies  were  performed  using 

the  same  potential  radius  8.5155  A  and  Duesbury  factor 
-1 

0.25a  as  for  small  clusters.  A  vacancy  "platelet"  was 
created  in  the  perfect  lattice  by  removing  the  desired 


number  of  adjacent  atoms  from  a  (111)  plane.  The  system 
was  then  relaxed  to  the  minimum  energy  configuration, 
which  took  from  20  to  100  iterations. 

Various  numbers  of  vacancies  were  removed  from  a 

(111)  plane  in  the  shapes  shown  in  Figure  7-2.  The  number 

of  vacancies  in  the  hexagonal  loop  in  a  (111)  plane  are 

given  by  3n(n+l)  +  1,  where  n  is  an  integer.  For  n  =  0 

one  gets  the  trivial  case  of  a  single  vacancy;  n  =  1  gives 

3 

a  platelet  with  seven  vacancies.  Since  Lam  ^  ^  has 
found  that  the  lowest  energy  configuration  of  a  nine 
vacancy  cluster  is  not  a  loop,  larger  sizes  were 
investigated.  The  two  smallest  hexagons  formed  by  ten  or 
more  vacancies  have  19  and  37  vacancies.  The  triangular 
shapes  are  given  by  n(n+l)/2;  sizes  of  10,  15,  21,  28  were 
calculated.  The  number  of  vacancies  in  rhomboids  can  by 
specified  by  m  x  n,  where  m  and  n  are  integers.  Rhomboids 
of  sizes  4x4,  4x5  and  5x5  were  simulated.  The  energies  of 
three  irregular  shapes  were  calculated;  a  19  vacancy 
hexagon  with  a  side  removed  to  give  a  16  vacancy  cluster; 
a  21  vacancy  triangle  with  two  corners  removed  to  form  a 
19  vacancy  cluster;  and  a  14  vacancy  cluster  formed  by 
removing  the  acute  corners  of  a  16  vacancy  rhomboid. 

The  first  goal  of  the  calculations  was  to  vacancy 
platelet  with  the  fewest  vacancies  which  would  collapse. 
The  criterion  for  collapse  is  that  atoms  above  and  below 
the  loop  move  to  almost  first  nearest-neighbor  distance 


from  each  other,  while  remaining  almost  planar.  Contrary 

3 

to  the  results  of  Lam  et  ^  ,  all  the  shapes  that  were 
tested  collapsed,  including  the  19  vacancy  hexagon.  The 
displacement  field  above  the  19  vacancy  loop  is  shown  in 
Figure  7-2,  which  shows  that  atoms  just  above  the  center 
of  the  loop  are  displaced  by  about  half  the  distance 
between  (111)  planes.  The  planes  below  the  loop  are  not 
shown  but  reflect  the  behavior  above  the  loop,  such  that 
the  opposing  atoms  are  at  first-neighbor  distances  from 
one  another.  The  loop  is  a  Frank  loop  since  there  is  no 
evidence  of  shear.  The  37-vacancy  loop  qualitatively 
collapsed  in  the  same  way. 

The  displacement  field  above  the  center  of  the  19  and 

37  vacancy  loops  is  plotted  quantitatively  in  Figures  7-4 

30 

and  7-5.  The  plot  is  taken  from  Ohr  and  shows  the 

predictions  of  anisotropic  elasticity  theory  for  the 
displacement  fields.  The  periodic  boundaries  used  in  the 
simulation  force  the  displacement  to  zero  midway  between 
images  of  the  loops,  although  this  zero  can  in  principle 
be  extended  to  large  radii  using  large  simulation  volumes. 
In  the  19  vacancy  loop,  the  displacement  field  appears 
have  converged  with  increasing  system  size  at  short 
distances  above  the  loop.  At  z/R  =  2.5,  the  displacement 
field  calculated  in  this  work  differs  by  at  least  40 
percent  from  anisotropic  theory  predictions.  The  37 
vacancy  results  are  inconclusive  here  because  the  field 
has  not  converged  with  a  system  size  of  N  = 


4000. 


Figure  2“1*  Displacement  in  z-direction  Above 
19  Vacancy  Loop 


The  lines  indicate  the  predictions  of 

30 

anisotropic  elasticity  by  Ohr  for  various 
metals.  The  symbols  are  the  results  of  this  work 
for  volumes  containing  the  specified  number  of 
atoms  before  vacancies  were  created.  R  is  the 
loop  radius;  b  is  the  Burgers  vector  a[lll]/3. 


Figure  7-5:  Displacement  in  z-directxon  Above 
37  Vacancy  Loop 


The  lines  indicate  the  predictions  of 

30 

anisotropic  elasticity  by  Ohr  for  various 
metals.  The  symbols  are  the  results  of  this  work 
for  volumes  containing  the  specified  number  of 
atoms  before  the  vacancies  were  created.  R  is  the 
loop  radius;  b  is  the  Burgers  vector  a[lll]/3. 
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The  triangular  vacancy  platelets  invariably  formed 
stacking-fault  tetrahedra  (SFT).  Figure  7-5  shows  three 
adjacent  (111)  planes,  including  the  plane  from  which  the 
vacancies  were  removed.  The  nearest  neighbor  atoms  one 
plane  above  the  circle  atoms  after  displacement  are 
connected  to  form  triangles.  The  triangles  which  contain 
the  circled  atoms  are  inverted  in  the  region  of  the 
removed  atoms  which  indicates  the  stacking  fault.  The 
displacement  of  the  atoms  after  relaxation  is  shown  more 
clearly  in  Figure  7-6,  which  is  a  cross-section  of  an  SFT 
created  by  21  vacancies.  The  stacking  fault  is  indicated 
by  the  S-labelled  triangles;  the  stair-rod  jog  by  the 
J-labelled  squares.  Note  that  the  atoms  below  the 
original  vacancies  have  displaced  only  slightly,  whereas 
the  atoms  above  the  vacancies  have  moved  about  4/5  of  the 
distance  between  (111)  planes.  In  this  view  the  atoms 
which  were  displaced  the  most  form  a  triangle  enclosed  by 
the  S-triangles  and  J-squares,  and  the  atoms  to  the  left 
and  right  of  the  triangle  are  only  slightly  displaced. 

The  collapsed  rhomboids  investigated  in  this  work 
resulted  in  asyraetric  displacement  fields  which  can  be 
thought  of  as  two  adjacent  SFT  which  are  inverted  in 
relation  to  one  another.  In  Figure  7-8  we  see  a  (IlO) 
plane  cross-section  of  a  25  vacancy  rhombus.  The  dashed 
line  indicates  the  position  of  the  original  vacancies; 
note  that  the  direction  of  the  collapse  is  opposite  on 


Three  planes  at  the  face  of  a  28-vacancy  SFT 
are  shown.  The  squares  are  the  atoms  in  the  plane 
from  which  the  vacancies  were  removed;  the  circles 
and  triangles  are  on  the  planes  below  and  above, 
respectively. 


Figure  7-7:  Cross-Section  of  21-Vacancy  SFT 


The  position  of  the  cross-section  is  shown  in 
the' caption.  The  solid  lines  connect  nearest 
neighbor  atoms  after  relaxation;  the  dashed  lines 
show  the  unrelaxed  positions  where  it 
significantly  differs  from  the  relaxed 
configuration.  The  S-triangles  and  J-squares 
indicated  the  stacking  fault  and  stair-rod  jog, 
respectively. 


i 
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either  side.  This  displacement  can  be  explained  by  noting 
that  the  25  vacancy  rhombus  can  be  constructed  from  two  15 
vacancy  triangles  (one  five  vacancy  edge  of  the  triangle 
is  counted  twice).  If  these  triangles  were  separated  but 
had  the  same  orientations  as  in  the  rhombus,  they  would 
form  two  SFT  with  opposite  orientations.  The  SFT  are  then 
merged  into  the  rhombus,  but  at  the  intersecting  edge  an 
intermediate  displacement  field  is  formed.  In  Figure  7-9 
three  adjacent  (li2)  planes  are  shown.  Notice  this 
displacement  field  is  qualitatively  identical  to  the  field 
around  a  similar  cross-section  in  the  19  vacancy  hexagonal 
loop  in  Figure  7-2.  This  is  not  unsurprising  when  it  is 
noted  that  the  25  vacancy  rhombus  can  be  formed  from  a  19 
vacancy  hexagon  if  three  vacancies  are  added  to  each  of 
two  opposite  edges  of  the  hexagon.  The  atoms  in  rhombus, 
however,  experience  a  shear  which  is  not  present  in  the 
hexagon.  In  Figure  7-10  three  adjacent  (111)  planes 
including  the  plane  in  which  the  vacancies  were  created 
are  shown.  The  dashed  triangles  are  drawn  between  the 
atoms  which  experience  the  large  displacement  on  either 
side  of  the  rhombus  as  shown  in  Figure  7-8.  The  atoms 
near  the  acute  corners  of  the  rhombus  are  displaced  in  a 
similar  manner  as  the  atoms  on  the  face  of  the  SFT  shown 
in  Figure  7-5,  but  there  is  evidence  of  shear  in  the 
rhombus.  In  the  center  of  the  rhombus  where  no  triangles 


are  drawn,  the  atoms  are  displaced  to  such  a  great  extent 


that  they  line  up  and  the  adjacent  planes  do  not  form  a 
stacking  fault,  but  are  in  intermediate  positions  and  are 
almost  unfaulted.  Hence,  although  the  rhombus  can  be 
considered  as  the  the  combination  of  two  SFT  or  as  a 
modified  hexagon,  we  see  that  the  final  structure  is 
different  than  the  structure  of  the  either  a  hexagonal 
loop  or  an  SFT. 

The  binding  energies  of  the  clusters  investigated  in 
this  work  are  given  in  Table  7-2  and  plotted  in  Figure 


7-11. 


Since  these  large  clusters  create  extensive 


displacement  fields,  it  was  necessary  to  use  large  systems 
to  achieve  or  at  least  indicate  whether  or  not  the 
calculated  binding  energies  have  converged  with  respect  to 
the  size  of  the  system.  The  polygonal  symbols  indicate 
the  binding  energies  in  an  864  atom  volume;  results  for 
larger  systems  are  shown  by  an  'x'  above  the  the  polygons. 
It  is  evident  from  Figure  7-11  that  the  SFT  are  the  most 
stable  over  all  the  cluster  sizes  and  configurations 
tested.  The  19-vacancy  "imperfect"  SFT,  formed  from  a 
21-vacancy  triangle  with  two  corners  missing,  has  a 
binding  energy  which  is  about  2  ev  greater  than  the 
binding  energy  of  the  "perfect"  19-vacancy  hexagonal  loop. 
The  rhomboids  are  also  more  stable  than  the  hexagonal 
loops,  but  less  stable  than  the  SFT. 
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Binding  Energy  [eV] 


Shape 

Cluster 

Size 

System 

Size 

Of 

Cluster 

Per 

Vacancy 

Hexagon 

16 

864 

5.79 

0.36 

19 

864 

6.81 

0.36 

19 

1372 

7.17 

0.38 

19 

2048 

7.41 

0.39 

19 

4000 

7.69 

0.40 

37 

864 

18.51 

0.50 

37 

2048 

20.53 

0.56 

37 

4000 

21.44 

0.58 

SFT 

10 

864 

3.68 

0.37 

15 

864 

7.27 

0.49 

19 

864 

9.34 

0.49 

21 

864 

11.86 

0.57 

21 

2048 

12.43 

0.59 

28 

864 

17.48 

0.62 

28 

2048 

18.39 

0.66 

28 

4000 

18.83 

0.67 

Rhomboid 

14 

864 

4.41 

0.35 

16 

864 

6.90 

0.43 

20 

864 

9.54 

0.48 

20 

2048 

10.18 

0.51 

20 

4000 

10.48 

0.52 

25 

864 

13.46 

0.54 

25 

2048 

14.23 

0.57 

25 

4000 

14.63 

0.59 

Table  7-2;  Binding  Energies  of  Large  Clusters 
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7.4  Discussion 


The  general  results  of  the  large  cluster  calculations 

in  this  work  follow  the  results  of  previous  calculations 

76,  77 

on  large  clusters  in  copper  .  It  is  found  that 

platelets  collapse  at  relatively  small  sizes,  starting 

with  as  few  as  10  to  15  vacancies.  The  triangular  vacancy 

platelets  form  stacking  fault  tetrahedra  (note  that  the 

stable  six  vacancy  cluster  for  copper  found  by  Lam  is  a 

3 

small  stacking  fault  tetrahedron  ). 

It  is  well  known  that  elasticity  theory  breaks  down 

near  a  dislocation  and  a  considerable  amount  of  simulation 

12 

work  has  gone  into  the  study  of  the  core  region  .  This 

is  probably  the  reason  for  the  discrepancy  in  the  results 

in  Figures  7-4  and  7-5.  Note  in  Figure  7-4  that  the 

difference  between  the  simulation  results  and  elasticity 

theory  is  fairly  large  when  compared  differences  in  the 

elasticity  theory  curves  for  different  metals.  Elasticity 

theory  depends  only  on  the  elastic  constants  of  a 

material,  and  the  elastic  constants  calculated  using 

Dagens'  potential  correspond  very  well  to  the  experimental 

15 

values  for  copper  .  Because  of  this  we  would  expect  the 
correspondence  between  the  the  copper  results  calculated 
from  the  simulation  and  from  elasticity  theory  to  be 
fairly  close  in  Figure  7-5,  if  the  elasticity  theory  were 
accurate  in  this  region.  This  is  especially  important 


since  anisotropic  elasticity  theory  is  used  to  interpret 

9,  11 

X-ray  diffuse  scattering  data 

Although  the  displacement  fields  around  large  loops 

have  been  previously  calculated  for  copper  using  atomistic 

76,  77 

simulation  with  empirical  potentials  ,  this  is  the 

first  calculation  of  the  binding  energies  of  these 

79 

clusters.  Bacon  and  Crocker  used  isotropic  elastictity 

to  calculate  the  energies  of  symmetrical  polygonal  loops 

and  found  a  circular  loop  to  be  the  most  stable,  but  did 

not  account  for  the  dislocation  core.  The  same  results 

were  found  when  the  theory  was  applied  to  nonregular 
80  81 

polygons  .  Brown  used  a  two-dimensional  analysis  of 

straight  dislocation  line  and  applied  it  to  calculating 

the  shape  of  loops,  and  concluded  that  the  variation  of 

core  energy  with  orientation  would  have  to  be  accounted 

for  to  find  that  hexagonal  loops  were  lower  in  energy  than 

circular  loops.  The  first  application  of  anisotropic 

82 

elasticity  theory  to  fee  loops  was  made  by  Bacon  et  al  , 

although  only  rhombus-shaped  loops  were  investigated. 

83 

Gaboriaud  and  Grilhe  calculated  the  relative  stability 
of  triangular  and  hexagonal  loops  and  SFT;  it  was  found 
that  SFT  were  more  stable  than  hexagonal  loops  of  about 
4400  vacancies  or  less  (this  corresponds  to  an  SFT  with  an 
240  A  edge).  These  calculations  accounted  for  the 
disassociation  of  the  partial  Frank  dislocation  around  the 


conclusions  of  Gabor laud  and  Grilhe  because  it  was  found 
that  the  SFT  were  the  most  energetically  stable  large 
clusters. 

Comparison  of  this  work  with  experimental  evidence  is 

8 

less  conclusive.  As  Kiritani  points  out  ,  the  relative 

binding  energy  of  clusters  is  not  the  the  only  factor 

which  determines  whether  or  not  a  cluster  will  be  formed. 

Unless  only  one  type  of  large  cluster  was  seen  under  all 

conditions,  it  would  be  difficult  to  determine  the  most 

energetically  stable  cluster  from  experimental  evidence. 

In  the  case  of  copper,  SFT  and  hexagonal  loops  have  both 
8,  84 

been  observed  .  The  sizes  of  the  collapsed  loops  also 

agree  with  X-ray  diffraction  studies  in  which  vacancy 

9 

clusters  as  small  as  about  10  A  in  diameter  are  found  , 
although  whether  the  cluster  was  a  loop  or  SFT  could  not 
be  determined. 

The  general  conclusion  of  this  work  is  that  loops  and 

SFT  can  form  in  copper  from  small  vacancy  clusters  with  as 

few  as  ten  vacancies.  It  is  interesting  to  note  that  the 

two  previous  simulations  of  these  large  clusters  in  copper 

yielded  essentially  the  same  conclusions,  even  though 

76,  77 

different  copper  potentials  were  used  .  This  is  one 

of  the  few  cases  where  all  the  potentials  for  the  same 
metal  have  given  similar  results  in  molecular  statics 


simulations. 


It  is  found  that  there  are  some 


discrepancies  between  the  results  of  this  work  and 

elasticity  theory,  and  it  is  concluded  that  the  latter  is 

in  error  for  clusters  of  this  size.  This  is  particularly 

30 

important  since  the  results  of  Ohr  are  used  to  interpret 

9,  11 

diffuse  X-ray  scattering  observations 
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Chapter  8 
Silver  Properties 


8.1  Introduction 


In  addition  to  the  copper  potential  used  in  the 

earlier  Chapters  of  this  work,  Dagens  also  developed 

13,  14,  15 

potentials  for  the  noble  metals  silver  and  gold 

Similar  to  the  copper  potential,  excellent  correspondence 

is  found  between  calculated  and  experimental  values  of 

15  17 

elastic  constants  and  the  phonon  spectra  for  silver. 
1,  3 

Lam  et  ^  calculated  the  formation  energies  of  single 

vacancies  and  interstitials,  binding  energies  of  some 

interstitial  and  vacancy  clusters,  and  the  migration 

1,  3 

energy  of  single  and  divacancies  in  silver.  Lam  ^  al 

calculated  the  migration  energy  the  single  vacancy  to  be 

0.57  eV,  which  within  twelve  percent  of  the  experimental 

4 

result  of  0.65  eV  .  A  large  discrepancy  was  found  in 
divacancy  migration  energies,  however,  because  the 
calculated  and  experimental  values  of  the  divacancy 
migration  energy  are  0.19  eV  and  0.57  eV,  respectively. 


A  review  of  the  annealing  literature  which  supports 
the  experimental  value  of  the  divacancy  migration  energy, 
however,  indicates  that  it  is  possible  to  question  the 


assignment.  The  divacancy  values  were  originally  reported 

in  two  separate  annealing  studies  of  quenched  silver  in 
85,  86 

1962  ,  The  primary  finding  of  both  investigations 

was  that  an  annealing  stage  exists  at  '■250K  which  anneals 

out  about  sixty  percent  of  the  residual  resistivity  caused 

by  the  defects  and  which  has  an  effective  activation 

energy  of  0.57  ev.  At  the  time  it  was  believed  that 

monovacancies  in  silver  had  a  migration  energy  of  about 

0.82  ev,  so  the  0.57  ev  value  was  assigned  to  divacancies. 

(The  single  vacancy  migration  energy  was  later  found  to  be 
87,  88 

about  0.65  ev  ).  The  second  major  finding  was  that 

the  annealing  followed  second  order  kinetics  throughout 

sixty  to  ninety  percent  of  the  annealing  stage  at  250 

K.  Second  order  kinetics  imply  that  two  of  the  same  type 

of  defect  are  involved  in  a  reaction.  It  was  suggested 

that  divacancies  were  combining  to  form  tetravacancies 

86 

through  most  of  the  annealing  in  the  stage  .  Following 

the  second  order  kinetics  was  a  faster  annealing  process, 

i.e.,  one  that  approached  first  order,  which  could  be 

explained  by  the  divacancies  going  to  clusters  which  had 

already  been  formed.  The  general  picture  which  emerged 

was  that  the  annealing  was  dominated  by  divacancies 

combining  with  each  other  during  most  of  the  annealing 

stage,  followed  by  divacancies  finally  combining  with 

clusters.  Tetravacancies  are  still  present  to  some  degree 

86 

at  the  end  of  the  annealing  stage 
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Because  of  the  resistivity  remaining  at  the  end  of 

the  stage,  it  is  possible  that  these  conclusions  may  be 

contradictory.  At  the  end  of  the  annealing  stage  at  250K, 

approximately  forty  percent  of  the  original  residual 

resistivity  remains.  If  most  of  original  defects  are  now 

members  of  tetravacancies,  the  resistivity  per  vacancy  of 

a  tetravacancy  would  be  about  forty  percent  of  that  of  a 

single  vacancy.  It  is  generally  assumed,  however,  that 

the  resistivity  per  vacancy  is  independent  of  the  cluster 

size  at  least  for  small  clusters  which  has  been  verified 

71,  89 

in  experiments  on  some  metals  .  Another  possible 

explanation  is  that  at  the  end  of  the  stage  the  defects 
are  members  of  larger  clusters  which  are  more  likely  to 
have  a  smaller  resistivity  per  vacancy.  At  the  onset  of 
first  order  kinetics  when  divacancies  start  combining  with 
clusters  instead  of  other  divacancies,  however,  the 
vacancies  must  be  partitioned  between  divacancies  and 
tetravacancies  since  larger  clusters  are  just  starting  to 
form.  At  this  point  at  least  thirty  percent  of  the 
original  resistivity  has  already  annealed  out;  iepending 
on  the  original  vacancy  conditions,  the  tetravacancies 
once  again  must  have  significantly  less  resistivity 
contribution  per  vacancy  than  single  and  divacancies  do. 

Another  interesting  result  reported  by  Antesburger  et 

al  also  leads  one  to  suspect  that  the  accepted  assignments 

87 

may  not  be  complete  .  Antesburger  et  a^  performed 


*• 
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annealing  studies  of  pure  and  doped  electron-irradiated 
silver.  The  purpose  of  doping  is  to  create  large  numbers 
of  vacancy  sinks  so  that  the  vacancies  will  interact 
mainly  at  the  sinks  instead  of  with  one  another.  This 
process  prevents  mobile  vacancy  clusters  from  forming 
because  the  vacancies  are  annihilated  at  sinks  and 
interstitial  clusters  or  absorbed  by  large  immobile 
vacancy  clusters  before  they  can  encounter  another 
vacancy.  Since  electron- irradiat ion  experiments  usually 
only  produce  single  vacancies  and  interstitials,  vacancy 
clusters  should  not  influence  the  calculated  activation 
energies.  It  was  found  that  the  activation  energy  for  the 
annealing  of  the  doped  silver  was  constant  throughout  the 
annealing  stage,  as  expected.  In  undoped  silver,  however, 
the  results  were  difficult  to  explain.  At  high  initial 
defect  densities,  behavior  similar  to  the  doped  silver  was 
found,  i.e.,  the  activation  energy  was  constant.  At  low 
initial  defect  densities,  the  activation  energy  was  found 
to  vary  from  a  maximum  of  0.65  ev  to  a  minimum  of  about 
0.57  ev.  Since  divacancies  are  created  by  a  second-order 
reaction  (the  combination  of  two  single  vacancies),  this 
is  surprising  because  one  would  expect  divacancy  effects 
to  increase  with  defect  density. 

The  general  procedure  in  annealing  studies  is  to 
attempt  to  explain  annealing  behavior  using  the  simplest 
possible  model.  If  a  single  defect  (i.e.,  single  vacancy) 


cannot  explain  the  results,  divacancies  and  occasionally 

trivacancies  are  introduced.  This  was  found  to  be 

4 

necessary  in  gold  and  aluminum  .  Often  there  is  support 

from  other  types  of  experiments  such  as  self-diffusion  and 

positron  annihilation  investigations  which  provide 

additional  information,  but  even  in  these  cases 

interpretations  can  be  incorrect.  Gold,  for  example,  is 

the  most  extensively  investigated  metal  using  annealing 

studies  with  investigations  dating  from  the  1950's.  The 

conclusions  of  an  extensive  fitting  study  performed  on  the 

90  4 

available  data  in  1978  were  generally  accepted  ,  and  yet 

91,  68 

later  found  to  be  incorrect 

Because  of  these  questions,  we  felt  it  would  be 

interesting  to  investigate  the  silver  potential  of 
15 

Dagens  further.  This  involves  the  identification  of  the 
migration  behavior  of  the  tri-  and  tetravacanc ies  and 
primarily  the  calculation  of  migration  energies.  Compared 
the  copper  work  in  Chapters  5,  6,  and  7,  fewer  atoms  were 
used  in  the  silver  study  for  both  the  small  defect  and 
thermodynamic  property  calculations,  and  no  large  defect 
calculations  were  performed.  Nevertheless,  these  results 
provide  important  groundwork  for  future  work  and  allow  us 
to  make  tentative  conclusions  concerning  the  silver 
potential . 
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8.2  Thermodynamic  Properties 

As  discussed  in  Section  5.2,  it  is  not  possible  to 

calculate  the  pressure  during  the  simulation  since  we  do 

not  know  the  structure- independent  part  of  the  potential. 

Instead  of  calculating  silver  properties  at  three 

different  volumes  as  was  done  with  copper,  however,  a 

reference  pressure  P  was  calculated  and  used  as  a 

ref 

correction  to  the  simulation  pressure  at  all  volumes. 

Since  the  silver  interatomic  potential  was  derived  for  the 

lattice  constant  a  =  4.087A,  the  reference  pressure  was 

calculated  at  T  =  298K,  since  real  silver  at  this 

temperature  and  ambient  pressure  has  the  desired  lattice 

constant,  was  found  to  be  P  =20.8  GPa. 

ref 

The  results  of  the  thermodynamic  calculations  are 
plotted  in  Figure  8-1.  Thes  calculations  were  done  with 
N  =  108  atoms  and  a  cutoff  radius  of  1.5a  for  the 
potential.  In  general  the  data  is  well-behaved  except 
near  the  melting  point,  which  is  to  be  expected.  Table 
8-1  contains  calculated  thermodynamic  properties  and 
experimental  values. 

The  results  of  bulk  modulus  calculations  using 
molecular  statics  are  shown  in  Figure  8-2,  and  the 
calculated  bulk  modulus  is  given  in  Table  8-1.  The  bulk 
modulus  was  calculated  at  T  =  OK  by  setting  up  the  system 
with  a  given  volume  and  atoms  at  the  perfect  lattice  sites 


This  Work 


92 


Exp 


c  [kJ/kg-K] 

P 

251. 

235@298K,  269@820K 

T  [K] 
m 

-1300 

1235. 

[kJ/Kg] 

m 

107. 

104.2 

aV 

f 

10 

0.014 

0.05 

B  [10  Pa] 

-6  -1 

10.1 

10.3 

^[10  K  ] 

10.7 

19.5 

Table  8-1;  Comparison  of  Silver  Thermodynamic  Properties 

and  then  calculating  the  pressure.  These  calculations 

were  done  with  three  different  systems  which  contained  a 

specified  number  of  atoms  and  potential  cutoff  radius. 

The  discontinuity  shown  in  the  (N  =  256,  r  =  1.866a)  case 

C 

is  caused  by  the  atoms  crossing  the  cutoff  radius  at  that 
particular  volume. 


8.3  Small  Cluster  Migration 


In  the  small  cluster  calculations,  and  cutoff  radius 

of  1.5a  was  used  with  a  system  of  108  atoms.  Also,  the 

45 

Duesbury  factor  was  not  employed  in  the  silver 

calculations.  Because  of  this  we  first  wanted  to  compare 

1,  3 

the  published  values  of  defect  results  of  Lam  et  al  to 


ours  to  insure  that  there  was  reasonable  agreement.  The 
migration  and  binding  energies  of  the  most  stable 
configurations  of  some  of  the  vacancy  clusters  are  shown 
in  Table  8-2  for  systems  of  108  and  256  atoms  (before  the 
atoms  were  removed  to  create  vacancies).  Table  8-2 
indicates  that  there  is  fairly  good  correspondence  between 
the  values.  This  is  not  unexpected  since  the  defect 
properties  are  mainly  determined  by  the  short  range 
portion  of  the  potential.  We  also  confirmed  that  the 
configurations  listed  were  the  most  stable. 


This 

Work 

a 

Ene 

£21 

Lam  et  al 

N'=108 

N’=256 

1 

m 

E 

IV 

0.54 

0.55 

b 

E 

2V 

0.04 

0.05 

1 

m 

E 

2V 

0.19 

0.22 

b 

E 

3V 

0.47 

0.39 

0.44 

b 

E 

4V 

0.86 

0.82 

Table  8-2: 

Comparison  of 

Silver  Results 

1,  3 

i  - 

Lam 

et  al 

To  determine  mechanisms  of  migration  and  breakup  of 
clusters,  it  is  necessary  to  know  the  binding  energies  of 
the  possible  cluster  configurations.  In  Figure  8-3  five 


trivacancy  configurations  are  shown.  Configuration  3E  is 
not  a  trivacancy  in  the  sense  that  all  the  vacancies  are 
not  connected.  The  binding  energies  of  the  configurations 
are  given  in  Table  8-3.  Configuration  3A  is  definitely 
the  most  stable  having  a  binding  energy  of  about  0.25  ev 
larger  than  the  others.  Configuration  3C,  which  from  a 
migration  standpoint  is  between  B  and  D,  is  less  stable 
than  3B  and  3D.  Even  though  configuration  3E  is  not  a 
trivacancy  by  our  definition,  it's  binding  energy  is 
fairly  close  to  that  of  3C. 

The  migration  energies  between  some  of  the 
configurations  are  given  in  Table  8-4.  The  trivacancy 
migrates  by  transforming  from  configuration  3A  to  3B  and 
then  back  to  3A  at  a  different  site.  As  shown  in  Figure 
8-4,  the  transformation  is  accomplished  by  the  movement  of 
two  atoms--the  central  atom  of  the  3A  configuration  and 
one  adjacent  atom.  The  migration  energy  for  the  motion  is 
0.47  ev.  Defining  a  reaction  coordinate  _  such  that  it  is 
1.0  near  in  configuration  3A  and  -1.0  when  near  3B,  we  can 
plot  the  energy  barrier  as  shown  in  the  left  half  of 
Figure  8-5. 

The  breakup  of  the  trivacancy  can  occur  by  a 
transformation  from  configuration  3B  to  3E,  as  shown  in 
Figure  8-4.  The  reaction  coordinate  is  defined  here  to  be 
1.0  at  3B  and  -1.0  at  3E.  As  shown  in  Figure  8-5,  the 
barrier  height  between  3B  and  3E  is  about  0.28  ev  above 


silver  Trivacancy  Configurations 


Table  8-3:  Binding  Energies  of 
Silver  Vacancy  Clusters 


Migration  Energy  (ev) 


Transition 

N‘=108 

N'*256 

3A->3B 

0.47 

0.47 

3B->3E 

0.28 

4A->4A' 

0.38 

4A->4C 

0.48 

4A->4D 

0.52 

Table  8-4: 

Silver  Migration 

Energies 

Silver  Migration  Mechanisms 


the  energy  in  configuration  3B.  We  note  that  the  energy 

at  configuration  3E  is  still  somewhat  greater  than  that  of 

an  isolated  divacancy,  so  the  single  and  divacancy  are 

still  attracted  to  one  another  at  this  point. 

The  situation  for  the  tetravacancy  is  more 

complicated  than  that  for  the  trivacancy.  We  investigated 

the  four  configurations  shown  in  Figure  8-6.  As  shown  in 

Table  8-3,  configuration  4A  is  the  lowest  energy 

3 

configuration,  as  reported  by  Lam  e^  al  .  The  other 
configurations  have  much  lower  binding  energies.  Note 
that  4D  is  not  really  a  tetravacancy  since  all  the 

vacancies  are  not  connected.  Some  of  the  migration 
mechanisms  for  the  tetravacancy  are  shown  in  Figure  8-4, 
and  the  migration  energies  are  given  in  Table  8-4.  It  is 
interesting  to  note  that  the  stable  configuration  4A  can 
migrate  without  passing  through  a  higher  energy  metastable 
conf igurat ion. 

8.4  Discussion 

The  thermodynamics  properties  that  were  calculated  in 

the  simulation  agree  are  in  good  agreement  with  the 

experimental  results.  The  two  quantities  that 

substantially  differ  are  the  thermal  linear  expansion 

coefficient  ^  and  the  volume  change  on  freezing  .  A 

f 

possible  explanation  for  the  large  errors  in  these 

quantities  is  that  the  reference  pressure  P 

ref 


,  which 


accounted  for  the  cohesive  energy  of  the  electrons,  is  not 
realistically  a  constant  over  the  wide  temperature  range 


(and  phase  change)  we  simulated.  Errors  in  P  would 

ref 

significantly  affect  and  because  they  depend  of 

volume  changes. 

The  defect  calculations  in  this  work  indicate 

interesting  defect  behavior.  The  most  interesting  result 

is  that  the  tetravacancy  is  apparently  more  mobile  than 

the  trivacancy.  This  is  caused  by  the  fact  that  the 

tetravacancy  can  migrate  without  passing  through  a  higher 

energy  metastable  configuration,  while  the  trivacancy 

passes  through  3B  during  migration. 

By  combining  the  results  of  this  work  with  the 

1,  3 

previous  studies  by  Lam  et  ^  ,  we  can  begin  to 

understand  the  annealing  behavior  of  the  computer  silver 

we  are  simulating.  In  Table  8-5,  the  basic  defect 

properties  affecting  annealing  are  shown.  In  cases  where 

Lam  ^  a^  had  previously  calculated  values,  we  use  their 

dis 

numbers  which  are  more  accurate.  The  column  E  is  the 
disassociat ion  or  breakup  energy  of  the  cluster.  This 
value  is  the  difference  in  binding  energies  before  and 
after  the  breakup,  added  to  the  migration  energy  of  the 
most  mobile  constituent. 

We  can  infer  some  conclusions  from  Table  8-5.  First 
of  all,  the  divacancies  are  so  mobile  compared  to  the 
other  clusters  that  they  will  be  annihilated  at  sinks  or 


Defect 


(ev) 

E  (ev) 

1  i 

* 

- 

0.55 

- 

* 

* 

0,04 

0.19 

0.59 

* 

0.47 

0.47 

0.62 

* 

0.86 

0.38 

0.86 

Table  8-5:  Defect  Annealing  Values  in  Silver 


1,  3 

{*  from  Lam  et  al  ) 


combine  with  other  vacancies  almost  instantaneously  upon 
formation.  This  can  shown  by  estimating  the  jump  rate  J 

2 

of  the  divacancies.  The  jump  rate  is  given  by 


J  =  F  exp(-E  /kT) 

2  2V 

13  -1 

where  F  is  a  frequency  factor.  For  F  =  10  sec  and  T  = 

8  -1 

200K,  we  have  J  10  sec  .  With  a  jump  rate  this  large 
2 

the  divacancies  will  encounter  sinks  or  other  defects  very 
quickly.  Because  of  this,  the  divacancy  concentration 
will  be  very  low  and  their  direct  effect  on  the  calculated 
activation  energy  will  be  small.  Even  though  the 
divacancy  binding  energy  is  very  small,  they  will  never 
disassociate  before  being  annihilated  because 


disassoc lat ion  requires  single  vacancy  migration.  The 
indirect  effect  of  the  mobile  divacancy  is  that  the 
trivacancy  will  breakup  or  disassociate  before  interacting 
with  other  defects  because  its  disassociation  energy  is 
low. 

The  consequence  of  low  divacancy  and  trivacancy 
concentrations  is  that  tetravacancy  concentrations  will 
also  be  small  because  the  divacancies  and  trivacancies  are 
required  in  reactions  that  create  tetravacancies.  Under 
these  circumstances  second-order  reactions  (divacancy- 
divacancy,  trivacancy-trivacancy,  etc.)  between  clusters 
will  be  minimal.  At  the  beginnning  of  annealing,  some 
pentavacancies  will  eventually  form  and  subsequently  most 
of  the  annealing  will  be  the  absorption  of  single  and 
divacancies  by  these  clusters.  We  conclude  that  this 
annealing  scheme  does  not  explain  the  experimental  results 
on  the  annealing  of  quenched  silver  any  better  than  the 
previous  model,  since  the  reaction  under  these 
circumstances  will  be  first  order,  i.e.,  single  vacancies 
interacting  with  sinks  and  the  newly  formed  vacancy 
clusters. 

The  major  difficulty  with  the  annealing  model 
suggested  by  the  computer  simulation  results  is  the  very 
high  mobility  of  the  divacancies.  The  problem  is  not  that 
the  divacancy  concentration  is  so  low,  but  that  the 
trivacancies  disassociate  so  easily.  If  trivacancies  were 


relatively  stable  during  annealing,  they  would  be  able  to 

combine  with  each  other  in  a  second-order  reaction 

resulting  in  the  formation  for  hexavacancies.  The 

problems  we  pointed  out  in  the  annealing  model  suggested 

86,  85 

by  the  quenching  studies  were  based  on 

tetravacancies  having  a  much  smaller  resist  it ivity  per 
vacancy  than  single  and  divacancies.  It  is  likely  that 
hexavacancies  have  a  smaller  resistivity  per  vacancy  than 
tetravacancies,  and  hence  it  would  somewhat  more 
acceptable  if  the  second  order  reaction  involved  tri- 
instead  of  divacancies. 

The  results  of  these  silver  calculations  support  the 
general  conclusions  of  the  copper  work.  Notably,  the 
silver  potential  reproduces  the  thermodynamic  properties 
of  real  silver  better  than  the  copper  potential  did  for 
copper,  and  even  does  well  in  estimating  the  melting  point 
(although  the  density  at  the  melting  point  is  incorrect). 
The  low  temperature  results  are  very  good,  as  was  found 
with  the  copper  potential.  As  before,  we  find  that  the 
migration  behavior  of  the  silver  tri-  and  tetravacancies 
is  not  easily  predictable.  In  silver  it  is  found  that  the 
divacancy  is  more  mobile  than  single  vacancies  and  tri- 
and  tetravacancies,  but  the  tetravacancies  are  more  mobile 


than  trivacancies. 


Chapter  9 
Conclusions 

The  purpose  of  this  thesis  is  to  apply  atomistic 
simulation  techniques  to  the  study  of  defect  clusters  in 
pure  metals.  Specifically,  the  method  of  molecular 
statics  has  been  used  to  investigate  vacancy  clusters 
using  an  interatomic  potential  derived  from  first 
principles  for  copper.  The  migration  energies  of  small 
clusters  (four  vacancies  or  less)  have  been  determined  by 
considering  all  possible  paths  and  associated  barriers  for 
moving  the  lowest  energy  configuration  from  one  location 
to  another.  The  stable  configurations  of  large  clusters 
of  10-40  vacancies  have  been  studied  by  forming  a  vacancy 
platelet  of  a  given  size  in  a  (111)  close  packed  plane, 
allowing  the  entire  system  to  relax  to  the  lowest  energy 
configuration,  examining  the  resulting  displacement  field, 
and  checking  for  the  collapse  of  adjacent  atomic  layers. 

In  order  to  examine  the  applicability  of  Dagens' 
potential  for  calculating  other  properties  of  copper,  the 
method  of  molecular  dynamics  has  been  used  to  study  bulk 
properties  over  a  temperature  range  up  to  melting.  A 
parallel  calculation  was  also  carried  out  using  and 
empirical  potential  fitted  to  the  bulk  modulus,  the 


-2 


lattice  constant  at  OK,  and  point  defect  properties. 

These  results  provide  additional  evidence  for  the 

usefulness  of  these  potentials. 

The  major  findings  of  the  vacancy  migration  energy 

study  are  the  migration  energies  for  tri-  and 

tetravacancies,  namely,  0.56  eV  and  0.38  eV.  These  are 

appreciably  lower  than  the  single  vacancy  migration  energy 

of  ~0.74  eV,  which  is  reasonably  well  established  from 

4 

annealing  experiments  .  The  divacancy  energy  has  not  been 

experimentally  measured,  although  the  divacancy  is 

5 

believed  to  be  more  mobile  than  the  single  vacancy  . 
Essentially  nothing  definitive  can  said  about  the  tri-  and 
tetravacancies  in  copper  from  experiment.  Molecular 
statics  calculations  based  on  Dagens'  potential  give  a 
value  of  0.82  eV  for  the  single  vacancy  migration  energy 
and  0.55  eV  for  the  divacancy.  The  present  results 
predict  that  trivacancies  should  be  just  as  mooile  as 
divacancies,  and  the  tetravacancies  even  more  mobile. 
This  prediction  needs  to  be  verified.  An  attempt  has  been 
made  to  model  the  existing  annealing  experiments  for 
copper  by  using  the  binding  and  migration  energies  of 
small  vacancy  clusters  calculated  using  Dagens'  potential. 
It  is  found  that  the  effective  activation  energy  is  still 
close  to  the  single  vacancy  migration  energy  for  a  variety 
of  sink  and  initial  single  vacancy  concentrations. 

In  this  thesis,  a  study  was  also  made  of  tri-  and 
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tetravacancy  migration  in  silver,  for  which  a  first 
principles  potential  is  also  available.  It  is  found  that 
the  tri-  and  tetravacancy  are  also  more  mobile  than  the 
single  vacancy,  but  in  this  case  they  are  both  less  mobile 

than  the  divacancy.  Combining  the  copper  and  silver 

results,  it  appears  that  di-,  tri-,  and  tetravacancy 
clusters  are  generally  more  mobile  than  the  single 

vacancy,  but  that  the  precise  relative  order  can  vary  from 
case  to  case. 

The  major  findings  of  the  study  of  stability  of  large 

vacancy  clusters  are  observations  of  collapse  of  platelets 

with  as  few  as  ten  vacancies,  collapse  hexagonal  platelets 

into  Frank  loops,  and  the  collapse  of  triangular  platelets 

into  stacking  fault  tetrahedra  (SFT).  Both  Frank  loops 

and  SFT  have  been  observed  in  copper  by  electron 

8,  84 

microscopy  ,  and  the  from  X-ray  diffuse  scattering  the 

smallest  Frank  loop  was  found  to  be  about  10  A  in 
9 

diameter  ,  approximately  the  size  of  the  19-vacancy 

displacement  loop.  From  the  simulation  results,  the 

displacement  field  above  and  below  the  Frank  loop  has  been 

calculated  and  compared  with  linear  anisotropic  elasticity 
30 

calculations  ;  close  to  the  loop,  the  two  results  differ 
which  in  not  surprising  in  view  of  the  fact  that  linear 
elasticity  theory  is  not  expected  to  be  valid  at  short 
distances.  This  effect  could  be  important  since  the 
elasticity  results  are  used  in  interpreting  X-ray  diffuse 
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scattering  data  .  The  binding  energies  of  the  Frank 
loops  and  SFT  have  been  calculated  and  the  results 
indicate  that  SFT  are  the  most  stable  large  vacancy 
cluster  configuration. 

The  bulk  property  study  shows  that  both  Dagens'  and 

the  empirical  potentials  given  reasonable  results  at  low 

temperatures.  This  is  perhaps  to  be  expected  in  the  case 

of  Dagens'  potential  since  it  is  already  known  that  the 

elastic  constants  and  phonon  spectra  are  calculated 
14 

correctly  .  As  elevated  temperatures  both  potentials  are 
less  satisfactory.  Although  Dagens'  potential  is  better 
theoretically  justified,  it  is  incomplete  without  the 
structure-independent  term  for  the  purposes  of 
thermodynamic  calculations. 

This  work  has  shown  that  it  is  possible  to  simulate 
relatively  large  vacancy  defect  clusters  using  atomistic 
simulation  methods.  It  would  be  interesting  to  repeat 
these  calculations  for  a  different  interatomic  potential 
which  has  a  higher  stacking  fault  energy  to  determine  if 
the  large  defects  are  significantly  different,  perhaps 
forming  unfaulted  loops  instead  of  faulted  loops  and  SFT. 
Furthermore,  the  dynamical  behavior  of  the  large  clusters 
can  be  investigated  using  molecular  dynamics  to  observe 
the  formation  of  stacking  fault  tetrahedra  and  unfaulting 
of  Frank  loops.  The  migration  energy  calculations  can  be 
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extended  by  performing  more  realistic  annealing  modeling 

93 

using  the  correlated  random  walk  method  ,  which  may  allow 
less  ambiguous  comparisons  to  be  made  with  experiment. 
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